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ABSTRACT 


Boost-phase intercept of a threat intercontinental ballistic missile (ICBM) is the 
first layer of a multi-layer defense. This thesis investigates the requirements and limita- 
tions of the U.S. space-based ICBM defense against North Korea, Iran and China by in- 
troducing an ICBM trajectory prediction, selecting an orbit for exo-atmospheric kill vehi- 
cles (EKV) and developing a hybrid guidance algorithm. The prediction of the ICBM 
trajectory takes the rotation of the earth and the atmospheric drag into account along with 
the gravitational forces and thrust. The threat ICBM locations, specifications and capa- 
bilities of the EKV and EKV carrier, and the capabilities of the space launch vehicle are 
analyzed to determine an appropriate orbit for the space-based intercept. The pursuit 
guidance, proportional navigation guidance and bang-bang guidance rules and their per- 
formances are investigated and simulated for three example ICBM threats in three- 
dimensional environment. The simulation results performances are compared and ana- 
lyzed for minimum miss distance, intercept time and total command effort. The guidance 
tules are combined to meet the mission requirements, resulting in a hybrid guidance algo- 
rithm, which uses different guidance rules for different stages of a boost-phase intercept 


scenario. 
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I. INTRODUCTION 


The National Missile Defense program was launched to protect the U.S. from bal- 
listic missile attacks. It is well known that the former Soviet Union has provided some 
countries with ballistic missiles and that some of these countries are upgrading or devel- 
oping their own missile systems [1]. The report of the “Commission to Assess the Ballis- 
tic Missile Threat to the United States” in 1998 showed the perimeter of the threat to 
Congress [2]. The report states that some of these countries have or are about to gain the 


ability to attack the U.S. with their missile systems. 


This study will investigate the requirements, limitations and the performance of 
the space-based, boost-phase intercontinental ballistic missile ICBM) defense. Boost- 
phase defense is investigated because previous studies showed that boost-phase defense 
has significant advantages when the electronic countermeasures are considered [3]. 
Space-based ICBM defense has been considered by the U.S. government since the begin- 
ning of the Strategic Defense Initiative [4]. It is desirable to have a space-based defense 
system to kill an ICBM in early stages of its trajectory [5]. The overall space-based 
ICBM defense scenario is illustrated in Figure 1. In the scenario, several EKV carriers are 
placed on a low-earth orbit to cover all possible ICBM threats. The ICBM launch is de- 
tected by the IR sensors in the geosynchronous orbit and tracked by RF sensors as shown 
in Figure 1. The track information is fused, and the command and control of the system is 
facilitated via the communication links between the assets as shown. After a launch deci- 


sion is made, the EKV is launched from the appropriate EKV carrier onto the orbit. 


EKV carriers 
IR Sensor orbiting the 
( Earth for Space- 
based defense 
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Figure 1. Overall space-based ICBM defense scenario. 


Ground-based defense requires that the exact location of the threat missile be 
known and that the interceptor be placed in a small pre-designated area [6]. The space- 
based defense can be designed to provide a continuous coverage over the globe with the 
selection of an appropriate orbit and by launching enough number of interceptors [7]. 
Deblois, et al. claim that the space-based interceptors have the same performance as the 
ground-based interceptors, and considering the cost-benefit ratio, they suggest that space- 
based ICBM defense is not desirable [8]. However, they make this comparison for fixed 
ICBM launch sites. In contemporary world, the ICBM threat may come from anywhere 
on the earth, and in order to cover the most possible threats we need a space-based de- 


fense system [4]. 


In this study, the exo-atmospheric kill vehicle (EKV) is considered to perform a 
hit-to-kill intercept. With the hit-to-kill intercept, the interceptor will smash into the tar- 
get at a high speed (approximately 7 km/s), which will create more destructive force than 
that created by the conventional explosives with the same mass [8]. In order to achieve 
this goal, an advanced and successful guidance algorithm is required. In this study, dif- 
ferent types of guidance rules are examined and a hybrid algorithm is developed for the 
EKV. 


A. THESIS OBJECTIVE 
The main objective of this thesis is to investigate the space-based, boost-phase 
ICBM interception. This includes the modeling of the ICBM motion, selecting an orbit to 


place the EKV carriers on and investigating the guidance laws to achieve the intercept. 


We developed a mathematical model to define the motion of an ICBM by taking 
the earth’s rotation and the atmospheric drag into account. Three different ICBM models 
having different mass fractions are generated for three different example launch points in 
North Korea, China and Iran. The initial launch parameters (the elevation launch angle 
and the azimuth angle) are estimated by using the Lambert guidance for impulsive mis- 


siles. 


The threat ICBM locations, specifications and capabilities of the EKV and EKV 
carrier, and the capabilities of the space launch vehicle are analyzed to determine an ap- 
propriate orbit for the space-based intercept. The inclination angle and the right-ascension 
angle of the orbit are determined to provide a uniform and continuous coverage over the 
three example countries. The altitude of the orbit is determined using the maximum range 


of the EKV and the capacity of the space launch vehicle. 


The pursuit guidance, proportional navigation guidance and bang-bang guidance 
methods are investigated for the EK V. These methods are implemented in the guidance 
unit of the EKV in a three-dimensional environment. The three guidance methods are 
simulated in SIMULINK®, and the results are presented. 

B. THESIS OUTLINE 

Chapter II of this study models the ICBM mathematically in the gravity field with 
the earth’s rotation and atmospheric drag. The mathematical model is implemented in 
MATLAB‘ and the launch parameters for a San Francisco attack from North Korea, 
China and Iran are delineated. In Chapter III, the required orbit for a successful space- 
based intercept against the ICBM launched from North Korea, China and Iran is derived. 
The orbital elements and the orbit requirements are introduced. Chapter IV introduces the 
principles of the pursuit guidance, proportional navigation guidance and the bang-bang 
guidance. The methodology of implementing these guidance rules in a three-dimensional 


model is explained. Chapter V presents the simulation results and compares the different 


guidance rules for the given scenarios. An analysis is based on simulation results con- 
ducted to identify the minimum intercept time with minimum guidance effort and mini- 
mum miss distance. Through trial end error the minimum miss distance is determined to 
achieve the hit-to-kill intercept. Chapter VI provides concluding remarks. Appendix A 


shows the flowchart of the simulation. 


Il. ICBM DYNAMICS WITH A ROTATING EARTH 


The threat intercontinental ballistic missile ICBM) launched from a hostile loca- 
tion is described in this chapter. The target ICBM in this study is a solid propellant, three 
stage, boosting missile reaching speeds above 6 km/s at the end of its boost-phase. The 
trajectory of the ICBM is derived as a function of the thrust that is generated by the solid 
propellant, the gravitational effects, the atmospheric drag and the rotation of the earth. 
The ICBMs are assumed to be launched from North Korea, China or Iran targeting San 
Francisco, California. 

A. BOOSTING TARGET MODELING 

In this chapter, the drag forces acting on the missile and the angular velocity re- 
sulting from the earth’s rotation are considered, and a close form solution is generated as 
a function of these forces, along with the thrust and the weight. 

1. Standard Atmosphere Model 

The atmospheric drag has a significant effect on the ICBM when moving in the 
atmosphere. The drag force is directly proportional to the air density, along with the 
cross-sectional area, drag coefficient and the square of the velocity. The air density de- 
creases proportional to the altitude and is measured in kg/m? in this study. The air den- 
sity of the atmosphere is modeled by using an exponential approximation [9]. This ap- 
proximation divides the atmosphere into two layers. The air density in the lower atmos- 
phere (below 9144 m) and upper atmosphere (above 9144 m), respectively, are given by 
[9] 

peai2esie A<9144 


i (2.1) 
Py =1.75228763e 4 = A> 9144 


where A is the altitude in meters and pis the air density in kg/m? . The U.S. Standard 
Atmosphere and the exponential approximation results are shown together in Figure 2. 
As shown in Figure 2, the exponential approximation for air density provides close match 


to the U.S. Standard Atmosphere. 
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Figure 2. Exponential approximation for air density. 


2. Boosting ICBM Description and Requirements 

For modeling the ICBM, the U.S. Peacekeeper missile is used as an example in 
this study because it is well published in open literature. The missile is a three-stage 
solid-propellant boosting ICBM capable of carrying up to 2270 kg of payload. The mass 
fraction, which is the ratio of the propellant weight to the total weight of the vehicle, dic- 
tates the velocity that the missile can attain by using its fuel [9]. The mass fraction of the 
ICBMs modeled in this study are 83%, 84% and 85%, respectively, to represent the 
North Korean, Chinese and Iranian ICBMs. Increasing the amount of propellant fuel does 
not necessarily increase the speed and range of the missile. The efficiency of the fuel can 
be understood from the mass fraction. The higher the mass fraction, the greater the speed 
and the range of the missile [6]. Since the ranges of the missiles from three countries 
mentioned above are different, the mass fractions and fuel amounts necessary are also 
different. The stage total masses, the propellant masses and mass fractions of the ICBMs 


are given in Table 1. 


Table 1. ICBM data matrix. 









































Stage 1 | Stage 2 Stage 3 | Payload | Mass Fraction 
Total mass (kg) 49000 27670 771 2268 
ss 83% 
s g Propellant mass (kg) 41640 | 23520 | 6554 0 
Total mass (kg) 53000 29670 8711 2268 
g 84% 
5 Propellant mass (kg) 45640 25520 7554 0 
Total mass (kg) 54000 30670 9711 2268 
a 85% 
§ | Propellant mass (kg) 46640 26520 8554 0 








The ICBM has to travel intercontinental distances to achieve its strategic mission. 


Hence, it needs to reach near orbital speeds to ensure these distances. The required veloc- 











ity to reach a given distance with a given launch angle can be calculated by using the 


rocket equation [9] 


(2.2) 


req — 





where gm is the gravitational constant (3.986 x10" m’/s? ), gis angular distance to be 
traveled, r, is the initial distance of the ICBM from the center of the earth, y is the eleva- 
tion launch angle, and a is the radius of the earth (6371 km). Using (2.2), the required 
speeds needed to reach the continental U.S. from Iran, China, and North Korea are shown 
in Figure 3. The launch angle is taken as 45 degrees, which ensures the maximum range 


for the required speed [9]. 
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Figure 3. Required speed for a given distance. 


As shown in Figure 3, the required speed increases as the distance between the 
target and the launch site increases. The speed required to reach San Francisco, California 
from Kilju-kun, North Korea is 7.49 km/s, which is in the region of a low earth orbiting 
vehicle speed. The required speeds for Xining, China and Busher, Iran launch sites are 
8.12 km/s and 8.65 km/s, respectively. Equation (2.2) assumes an impulsive missile, 
which is not true for a boosting missile. The speed requirements for a boosting missile 
will be less than for an impulsive missile because there will be thrust acting on the mis- 
sile for several minutes. In addition, the rotation of the earth and the drag forces are not 
considered in this rocket equation. 

3. Boosting ICBM Mathematical Modeling 

In this section, we derive the mathematical model for a boosting ICBM that takes 
the earth’s rotation and the atmospheric drag into account. Kashiwagi derives a full 
mathematical model for re-entry vehicles, where the non-accelerating vehicle is released 
from space [10]. We adopt his derivation for ground-based boosting ICBMs by adding 
the thrust force generated by the solid propellant fuel. 


Before starting to develop a closed form solution for predicting the trajectory of 


the ICBM, we briefly define the coordinate systems. The common coordinate system 
8 


used in this derivation is the Earth-centered Earth-fixed (ECEF) coordinate system. The 
earth is assumed to be a perfect sphere in this coordinate system. The local coordinates, 
such as geodetic and topo-centric horizon, are transformed to the ECEF coordinate sys- 
tem. 

The ECEF coordinate system is an orthogonal Cartesian coordinate system that 
takes the center of the earth as the origin. The x-axis passes through Greenwich, y-axis 
passes through E90, and the z-axis passes through the North Pole. The angular velocity 


resulting from the earth’s rotation is added to the motion of the missile. 





North Pole 





Missile 
, > body 
centerline 
w90 
> 
y: 
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ny 
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Figure 4. ECEF and other coordinate systems. 


The topo-centric coordinate system assumes a locally flat earth. The Launch an- 
gles are commonly given in this coordinate system [6]. The rotation of the earth is de- 


noted by was shown in Figure 4. 


In the derivation of the differential equation for the motion, the ICBM is consid- 
ered a point mass. The force on the point mass, based on using Newton’s second law, is 


given by 

F=ma (2.3) 
where @ is the acceleration in m/s? and m is the point mass in kg. The unit of F is speci- 
fied in N (Newtons). 

The angular velocity vector of the earth, denoted by @, is given by [10] 

@=on (2.4) 
where fi is the unit vector in the z direction, and || = @ = 2m rad/day. In Figure 4, E isa 
point on the earth, and the vector S representing the ICBM missile body from point E to 
point P is given by 

S=xityj+zk (2.5) 
where x, y and z are the distance measures in ECEF coordinates, and 7, J and & are the 


unit vectors in the direction of x, y and z, respectively. The velocity of the ICBM with re- 


spect to point E can be defined as the time derivative of S : 


Va asi + jth (2.6) 


a5 


The acceleration of the point mass with respect to point £ can also be written as 


time derivative of the velocity: 
Pelee eee (2.7) 
dt 


The acceleration acting on the point mass is called the apparent acceleration and is de- 


noted by 


a, =a+a,+a, (2.8) 


P 
where dis the relative acceleration, d, is coincident point velocity, and a, is the coriolis 


acceleration [10]: 
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a, = @x(@x?F) (2.9) 
a, =20xV (2.10) 
By substituting (2.7), (2.9) and (2.10) into (2.8), we obtain the resultant apparent 
acceleration that is acting on the point mass [10]: 
a, = 31 + jj +2k + @x(@xF)+20x0 (2.11) 


The rotation of the earth @ can be defined as a function of ECEF coordinates, 
which will make it easier to solve for a closed form equation, and in doing so, we will 


keep a common coordinate system for calculations [10]: 
@=asin(u)k +@cos(u) 7 (2.12) 


The vector 7 , which defines the position of the ICBM, can also be written in the 
ECEF coordinate system as follows [10] 
F=xi+yj+zk (2.13) 


The cross products appearing in (2.11) are solved as shown below [10]: 


@x(@xF)=-0"xi +0" | (yoos(z )+zsin(u ))eos (ys )-y]i 
+0°[ (yeos(s1)+zsin(4e ))sin( L)-z zk (2.14) 
jn 


@xV = o(2cos(u)—psin(w))i +xesin(“) zacos(u)k 


where sz is the geodetic latitude of the ICBM. Substituting (2.14) into (2.11) yields the 
closed form of the acceleration acting on the ICBM [10]: 
= [z + 2o(2 cos(w)—jsin (#))]é 
+/+ 2kosin(u)+o° (ycos(u) +zsin(z2))cos(z) | (2.15) 
+[2-2kec0s(u)+ o (vcos(u) +zsin(w))sin(s2) |& 


Equation (2.15) concludes the derivation of the acceleration acting on the ICBM due to 


Newton’s second law (2.3). 


Now the left-hand side of (2.3), which contains the forces acting on the ICBM, is 
derived. The four major forces acting on the ICBM are the weight (W), the drag force 
(F,), the lift and the thrust (7). The ICBM examined here is a cruciform missile, hence 
the lift on this missile is very small and can be neglected. The total force on the ICBM is 


the sum of all remaining force vectors as given by 


F= 


om 


/+T+W (2.16) 





Figure 5. Forces acting on missile. 


The thrust vector 7 is along the velocity vector of the missile as shown in Figure 
5. The magnitude of the thrust vector is proportional to the propellant used and the spe- 


cific impulse of the missile /,,, given by 


== (2.17) 


where W represents change in mass of the propellant. In other words, it is the fuel spent 
ina second. The specific impulse is a measure of effectiveness of the fuel in units of sec- 
onds. The more capable missiles have a higher value of /,, [9]. The ICBM examined in 
this study is assumed to have a specific impulse of 300 s. The derivative can also be de- 
fined as the product of the gravitational acceleration and the change in mass [6] as fol- 


lows: 


W =mg (2.18) 


where g is the gravitational acceleration of the earth in m/s? . The gravitational accelera- 


tion of the earth varies inversely proportional to the square of the altitude of the missile 


g= GM (2.19) 


r 





where G is the gravitational constant of the earth, which is approximately equal to 


6.7154x10""' Nm/kg? [11], M is the mass of the earth, which is approximately equal to 
5.98x10™ kg [12], and r is the magnitude of the vector 7 as shown in Figure 4. 
When we substitute (2.18) and (2.19) into (2.17) and solve for the thrust, we have 


aclosed form solution of the thrust vector in the direction of the velocity vector V as 


given by 





(37 + i + 2h) (2.20) 


as > ey “a m I] 
7 pees [eae] 28h 
Ir I 
The drag force acting on the ICBM is a function of atmospheric density (p ), gravita- 
tional acceleration of the earth g, the ballistic coefficient 2 , and the velocity of the 


ICBM. The drag force is defined in the opposite direction of the velocity vector [9]: 


F, POP si 2) (2.21) 


The weight vector of the ICBM is towards the center of the earth and is defined as the 


product of the mass and the gravitational acceleration g as given by 


W=- 








most +) BOM sia se) (2.22) 
Now we can write the left-hand side of the equation, which defines the forces acting on 
the missile. Substituting (2.20), (2.21) and (2.22) into (2.16) gives the total forces acting 


on the ICBM as follows: 








rot noPl., a 
r 


|r| 28 r 








.GM, x mp\V|. GM _), 
«(0S a lm Bj (2.23) 











[ats aad 
” I7 2p r 








Substituting (2.23) and (2.15) into (2.3), and the rearranging the equation for the 
ECEF components, yields the differential equations of the acceleration as given by 
pI. tnGM 


¥=-2o(Zcosw—jsiny)+o° x 3 ps (xt 


pV |. GM GM 
2p. Pr? my” 
pl]. GM _ GM 

—7z+ I 


2B as r . (P| 9 














py (2.24) 








2exsin uz +o ysin? 4+ @"zsin "cos 











2 =2excos + @zcos’ “—w' ysin “cos 


The state vector of the ICBM is defined as a function of its position and velocity 


and denoted by 
X=[x y 2 V, V, Vale yp z ep ey (2.25) 
The evaluation of the change in the state of the missile in discrete time is given by 


(2.24). Representing these equations in ¥ = FX(t —h i format will require defining the 


state transition matrix F. The transition matrix is formed by using (2.24). The parameter ¢ 


is the time of interest and ¢, is the initial time [13]: 


























0 0 0 1 0 0 
0 0 0 0 1 0 
0 0 0 0 0 
mel, 
@- Gl 0 0 eo _ PB 2osinu  —2a@cosu 
F= r mV | 2p 
Bes : , mgl,, pg 
0 @ sin u- —@ sin L1.cos —2a@sin eats 
sin ft A iN {LCOS LL in LL mV 2B 
9 Ss mgl,, pg 
0 —@ sin cos @ cos(4)— 2a@cos 0 aH 
sin ££cos (2) 3 $f mlV| 2B 
(2.26) 
where 


@ : earth rotation rate (27/day) 

G : gravitational constant of earth (6.67 x10") 
M : mass of earth (5.98x10™ kg) 

r : distance of the target from earth's center 

m : fuel consumption of the target ICBM (kg/s) 
I,, : specific impulse of the target ICBM 

m : total mass of the target ICBM 

|| : magnitude of the velocity of the ICBM 

p . atmospheric density (kg/m*) 

B : ballistic coefficient, and 

i: geodetic latitude of the target ICBM 


The first order Markov model for the target’s state transition uses the state transi- 


tion matrix F. The differential equation derived in the preceding discussion is the mathe- 


matical model for the motion of the boosting ICBM in the gravity field. It is very cum- 


bersome to derive a closed form solution for this equation. Instead of deriving a closed 


form solution, we apply numerical integration using MATLAB". The numerical integra- 


tion method guarantees accurate results but demands high computational power. The mo- 


tion of the ICBM is simulated using the discrete-time implementation of the transition 


matrix and the state space vectors as follows 
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X(k+At)=X(k)+FX(k)At (2.27) 


where At is the time step. 

4, Initial Values of ICBM Launch Angles 

The azimuth launch angle C is measured in the topo-centric coordinate system 
from the North Pole to the tip of the missile. The measurement is taken from north to 
east. Consider the triangle ABC shown in Figure 6 as the launch geometry, where C is the 
launch point and B is the target location. The North Pole is denoted by A in this geome- 
try. The initial azimuth launch angle of the ICBM is denoted by C. 


North Pole 






Launch 


Point -:-> Target 





E90 


Greenwich 


Figure 6. Launch parameter geometry (After [15]). 


The lowercase letters in Figure 6 correspond to the angles subtended by the arcs 
measured form the center of the earth. The capital letters in Figure 6 correspond to the 
angles formed by the intersecting arcs. Note that angle C is the azimuth launch angle. Us- 


ing the law of sines, we can obtain the relationship among the angles as 


sina _ sine 








= 2.28 
sind sinC ( ) 
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By rearranging (2.28), the azimuth launch angle can be written as 


Bes en (2.29) 
sinacosC 


The denominator of (2.29) can be shown to be a function of angles a, b and c as given by 


sina og C= SO8e Ons cas, (2.30) 
sinb 


from the law of cosines: cosc = cosacosb+sinasinbcosC. By substituting (2.30) into 
(2.29), the azimuth launch angle is obtained as 


tanC= sincsinbsin A 231) 
cosc—cosacosbh 


The angles A, b and c are defined by the target location and the launch site location. The 
only unknown in (2.31) is the cosa term, which in turn can be written as a function of 


the known parameters by using the following law of cosines: 
cosa =cosbcosc+sinbsinccos A (2.32) 


The known angles mentioned above are defined using the geodetic locations of 


the target and the launch site as given by [15] 


A=A-A, 
a, a 

poke 2.33 
5 He (2.33) 

cane 


where 4 is target geodetic longitude, vis target geodetic latitude, 1, is launch site geo- 


detic longitude, and yz, is launch site geodetic latitude. 


Using (2.33), we can show the azimuth launch angle to be a function of geodetic 


locations and the cosa term as given by 


tan C= SOSH EOS Mi sin(A-A,) (234) 


sin 4¢—cosasin 1, 


The azimuth launch angles from Kilju-kun, North Korea, Xining, China and Bushehr, 
Iran to San Francisco (N37.76 W122) are calculated by using (2.34). The results are tabu- 


lated in Table 2 in which the azimuth launch angles are defined from the North Pole east- 














ward. 

Table 2. Azimuth launch angles for launch attitude 
Launch Site North Korea China Iran 
Location (geodetic) N41-E129 N36-E101 N29-ES1 








Azimuth launch angle 
40.37 31.52 353.24 
(degrees) 











The elevation angle for the ICBM is calculated using the Lambert guidance. 








Lambert guidance will put the ICBM on a collision triangle that is moving in a gravity 
field. We will solve the Lambert’s problem using a numerical method. In this solution, 
we assume a flat earth and use topo-centric coordinates. The elevation angle of the ICBM 
is denoted by y and measured from the earth’s surface to the tip of the missile. The cen- 
tral angular distance to be traveled is denoted by ¢ and measured from the earth’s center 
between the target location and the launch location [9]. The radius of the earth is denoted 
by R.. Note that (2.2) gives the required velocity of the ICBM for a given distance. Since 


the launch point and target location are both on the ground, 
r,=a=R (2.35) 


is a valid statement. Using (2.35) we can rewrite the required velocity equation in closed 


form [9] as 


(2.36) 





The time of flight for the ICBM can be calculated by using the formula for the el- 
liptical travel as given by [9] 





tan y(I1—cosg)+(1—A)sing 


l-cos@ + cos(y +) 
Acos’ y cosy (2.37) 


i tiie 
Veg ©08(7) e-a] 





2cosy 


Af (2/a)-1]" 


ee (2/4) 1 
cos y cot(¢/2)—siny 


where / is a constant and is given by [9] 
2 R, 
(2.38) 


qa 


gm 


The central angular distance ¢ can be calculated using the position vectors of the 


launch point and the target location: 
(2.39) 





where 7; and 7; are the launch point position vector and target location position vectors, 
respectively, in the ECEF coordinate system. 
By substituting (2.38) into (2.36) and solving for the elevation angle, we obtain 
the minimum and maximum possible elevation angles as follows [9]: 


sing—, ea 


Yorn = tan™ ( (Ico) 
(2.40) 


Sean sin $+ ,/2(1—cos¢) 
Fa = cod) 


We will use the method described in [9] to find the elevation angle corresponding 
to the desired flight time tr» which is calculated by using (2.37). The flight time is cal- 
culated iteratively by using elevation angles between the y,,, and 7,,,, (see (2.40)) in or- 


der to reach a value that is satisfactorily close to the desired flight time as follows [9] 


(7, = PaV(te -t,) 






Visita? (2.41) 


where n is the index of iteration. The elevation angle y is computed for North Korea, 


China and Iran are 54.22°, 55.98° and 58.92°, respectively. 


Note that the Lambert Solution assumes an impulsive missile moving in free 
space. However, the ICBM model created is a three-stage solid-propellant missile, hence 
the actual elevation angles for a boosting ICBM moving on a rotating earth and in the at- 
mosphere are different from the above values. The elevation angle for a boosting missile 
should be above 80° to overcome the gravity force and avoid hitting the ground [9]. To 
find the accurate elevation angles, the three-dimensional motion simulation is run and the 
results are reported in the following section. 

B. SIMULATION RESULTS FOR THE ICBM MODEL 

The preceding analysis is implemented numerically and simulated by using 
MATLAB‘. The booster fuel consumption is simulated by decreasing the amount for 
every time step. The time step is selected as 0.05 s as a trade-off between the accuracy of 
integration and the available computational power. The initial state vector of the ICBM is 
denoted by its launch location and launch attitude. The launch attitude is defined by the 
azimuth launch angle and elevation launch angle, which are given in topo-centric coordi- 
nates and transferred to the ECEF coordinate system. The initial attitude of the missile is 
represented in the initial velocity vector of the ICBM. The location of the missile is given 
in the geodetic coordinate system and it is also transferred to the ECEF coordinate sys- 
tem. The location of the missile is represented in the position portion of the state vector. 
Having the initial state vector and the transition matrix meets the requirements to run the 
simulation. The launch angles for the ICBMs launched from North Korea, China and Iran 


to reach San Francisco are given in Table 3. 
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Table 3. Initial launch parameters for ICBMs. 























North Korea China Iran 
Elevation launch angle in de- 
85.81 85.15 84.1 
grees 
Azimuth launch angle in de- 
39.6 19.8 342.1 
grees from North to East 
Location in geodetic coordi- 
N41 E129 N36 E101 N29 E51 
nates 

















The initial state vectors in the ECEF coordinate system for a San Francisco attack 


for running simulations are as follows 


-3025933.19 m -983476.74 m 3506700.44 m 
3736716.22 m 5059549.23 m 4330414.40 m 
4179752.70 m 3744779.84 m 3088722.09 m 
Xy korea (0) = 5X china (0) = 5X jan (0) = 
. -9.10 m/s -3.11 m/s 9.76 m/s 
9.82 m/s 13.32 m/s 11.15 m/s 
13.04 m/s 11.70 m/s 10.22 m/s 
(2.42) 


Given the initial state vectors, the state vectors at a future time instant t+ At can be com- 


puted using the iterative relationship [10] 

X(t+At)= X(t)+ FX (t)At (2.43) 
where At is the iteration time step and F is the state transition matrix. Equation (2.43) is 
iterated in time along with the decreasing amount of propellant fuel. The rotation of the 


earth, the drag force and the thrust are also integrated in the simulation as outlined in 


(2.26). 
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The geometry of the trajectories of the ICBMs launched from North Korea, China 
and Iran are illustrated in Figure 7. The trajectory of the Iranian ICBM passes through the 


North Pole, which explains the reason why the azimuth launch angle is towards the north. 





ncisco 


Position in z-direction in km 


-5000 


S000 


Position in x-direction in krn Position in y-direction in km 


Figure 7. Trajectories of the ICBMs launched from North Korea, China and Iran. 


The magnitudes of the velocities of the three ICBMs obtained from the simulation 
are shown in Figure 8. The velocity increases during the boost-phase, which lasts for 
about three minutes. After the boost-phase, the velocity starts to decrease because of the 
gravitational effect while the ICBM continues ascending. After the ICBM starts to turn 
downwards, the gravitational force acts in favor of the ICBM velocity and increases the 
velocity. This increase in the velocity persists until the ICBM enters the atmosphere. 
Once the ICBM enters the atmosphere, the drag forces become significant and decrease 


the ICBM velocity. 
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Figure 8. Magnitudes of the velocity vectors of the ICBMs throughout the entire flight. 
The earth’s rotation and the atmospheric drag are included 


The boost-phase, where the velocity increases very rapidly, is the most dynamic 
region for the ICBM. In taking a close look at the velocity graph during the boost-phase, 
Figure 9 illustrates the stage transitions of the ICBM. Remember that the required veloc- 
ity of the ICBM is proportional to the distance between the target and the launch point. 
Note that the maximum velocity at the end of the boost-phase has the increasing order in 
sequence of North Korea, China and Iran. This sequence is also true for the distance be- 


tween these launch points and San Francisco. 
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Magnitudes of the velocity vectors of the ICBMs during boost-phase. 


In addition, the required velocity is approximately 1 km/s less than that estimated 


using closed form rocket equation (2.2). The required velocities computed by both the 


rocket equation and the simulation are shown in Table 4. The simulation results are more 


reliable than the rocket equation results because they include the boosting ICBM condi- 


tions whereas the rocket equation results assumes an impulsive missile (initial velocities 


model). 


Table 4. 


Required velocities. 





Required Velocity (km/s) 




















Launch Site é - - Distance (km) 
Rocket Equation Results | Simulation Results 
North Korea 749 6.55 8700 
China 8.12 6.98 10700 
Tran 8.65 7.47 12500 

















The altitudes of the ICBMs are calculated by using the state vectors. Since the 


earth is modeled as a sphere, the altitude of the ICBM can be written as a function of the 


ICBM position: 





(2.44) 





where x, y and z are the position components of the ICBM in the | ECEF coordinate 


system, and R, is the radius of the earth. The altitudes of the ICBMs are shown in Figure 


10. An exponential 


crease is a result 0 
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increase in the altitude is observed during the boost-phase, which in- 


the increasing velocity. 
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Figure 10. Altitudes of the ICBMs during the entire flight. 


The total mass of the ICBM is shown in Figure 11. The total mass is the sum of 


the payload, the remaining propellant mass and the remaining canister mass. The total 


mass decreases as the ICBM burns the propellant fuel and releases the canister after the 


depletion of the propellant fuel. At the end of the boost-phase, only the payload is left 


having a mass of 2268 kg. 
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Figure 11. Total mass of the ICBMs as a function of the time from launch. 


After the depletion of the propellant mass in a stage, the empty canister is released 
and use of the propellant fuel in the next stage begins. The decrease of the propellant 
mass is shown in Figure 12. The propellant mass is used up in three minutes, which is the 


end of the boost-phase. 
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Figure 12. Propellant mass of the ICBMs during the boost-phase. 
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The overall simulation results are tabulated in Table 5. The burn out time, the alti- 
tude of the ICBM and the burn out velocities are tabulated for each stage transition. The 
maximum altitude is the utmost altitude of the ICBM during the entire flight. The total 


flight time denotes the time when the ICBM hits the target. 


Table 5. Simulation results for ICBMs from North Korea, China and Iran. 





Stage 1 Stage 2 Stage 3 





Launch Site 


Altitude (km) 
Velocity (km/s) 
Altitude (km) 
Velocity (km/s) 
Velocity (km/s) 


Altitude (km) 
Maximum Altitude (km) 


Total Flight Time (minute) 











North Korea | 28.73 | 1.44 | 119.34 | 3.9 281.66 | 6.55 | 1959 | 37.96 
China 28.24 | 1.50 | 112.94 | 4.02 | 262.2 6.98 | 2047 | 42.32 
Iran 24.31 | 1.52 | 85.56 4.1 182.02 | 7.47 | 1202 | 37.47 















































. SUMMARY 

In this chapter, a mathematical model for the ICBM motion is derived taking the 
earth’s rotation and the atmospheric drag into account. The three example ICBMs 
launched from North Korea, China and Iran have different mass fractions in order to 
reach San Francisco from different ranges. The entire flight of these example ICBMs are 


simulated by using the developed mathematical model. 
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II. ORBIT SELECTION FOR EKV CARRIERS 


In this chapter, we will describe the configuration and specifications of the space- 
based exo-atmospheric kill vehicle (EKV) carriers and the EKV itself, and determine an 
orbit in which to place these EKV carriers in order to intercept ICBMs launched from 
three example launch sites in North Korea, China and Iran prior to the release of the reen- 
try vehicles (RVs). The trajectories of the ICBMs were derived and simulated in the pre- 
vious chapter and the resulting trajectories are used to determine the orbit that is required 
for a successful intercept. We will show that for the given launch locations, we need a 
circular orbit with an altitude of 1000 km, an inclination angle of 43.5° and a right- 
ascension angle of 15.3°. 

A. INTERCEPTOR MISSILE MODELING 

In this study, we will model an EKV that conducts a hit-to-kill intercept. Hit-to- 
kill interception is selected because the damage applied in space is significantly more 
than the damage applied by conventional explosive interceptors. The EKV must hit the 


payload section of the ICBM in order to assure the desired hard kill [16]. 


Raytheon has developed a ground-based EKV and tested it successfully. The 
model of the EKV considered in this study will be based on the specifications of the Ray- 
theon EKV. We will use this model as a space-based EK V instead of a ground-based 
EKV [17]. The Raytheon’s EKV is shown in Figure 13. 


The IR detectors of the EKV are capable of determining the spatial location of the 
target and discriminating the target from the decoys [17]. The scope of this study is fo- 
cused on the boost-phase intercept in which it is assumed that the ICBM cannot deploy 
its decoys. However, this discrimination capability gives the EKV the flexibility to be 
used in mid-course as well. The EKV kills its target through the force created by the im- 
pact. 

The EKV has four thrusters for applying the guidance command generated by the 
guidance unit, as shown in Figure 13. In addition, the EK V has two attitude controllers 
which will cooperate with the divert thrusters and control the attitude of the EKV. Note 


that the EKV would tumble in outer space without these attitude controllers. The EKV 
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has onboard sensor optics to track the target and a guidance unit which improves the tar- 
get data refresh rate and decreases the delay between the tracking and guidance applica- 


tion. The optic sensors are protected by a sunshade cover, as shown in Figure 13. 
An ICBM has two different IR sources, including the plume and the body friction. 
The plume has a mid-wave (3—5 pm) infrared signature generated by high temperature 


emissions from the nozzle. The body friction produces a long-wave infrared signature 


(8-12 pm) [18]. The EKV has to have two different infrared sensor sets to track both of 


these sources. 


The EKV weighs 64 kg with a length of 139.7 cm and a diameter of 61 cm [17]. 
We added a 136 kg of booster to the EKV to give it the initial velocity when launched 
from the EKV carrier. The solid propellant in the booster is 100 kg, which makes the 
EKV’s total mass 200 kg. 
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Figure 13. Raytheon’s exo-atmospheric kill vehicle (From [17]). 


The EKVs are stored in a space-based EKV carrier, which is assumed to hold 10 
EKVs. The EKV carriers travel on an orbit that provides enough kinetic energy to allow 
the EKVs to succeed in destroying the ICBM before it delivers the reentry vehicles 
(RVs). Accomplishing a boost-phase intercept will decrease the amount of targets that 


need to be dealt with. 


The EKV is modeled using the same methodology that was introduced in Chapter 


II. The EKV boosts for 10 s and then moves into the gravity field with the guidance com- 
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mand. The guidance command is applied by the onboard divert thrusters in pitch and yaw 
axis. The guidance method and application will be discussed thoroughly in Chapter IV. 
B. ORBITAL ELEMENTS 

The laws that describe the motion of the planets in the solar system also describe 
the motion of a spacecraft orbiting around a planet. Johannes Kepler (1571-1630) derived 
three empirical laws to describe the motion of planets by using his observations. The mo- 
tion of any two bodies in space can be characterized by his laws [12]. 

1. Keplerian Orbit 

We introduce the Keplerian orbit and orbit perturbations in the following discus- 
sion. The shape, size and orientation of the orbit is defined by its orbital elements, which 


are illustrated in Figure 14 [19] and are described below briefly. 


Perigee is a point on the orbit that is closest to the earth’s center. Apogee is a 
point on the orbit that is farthest from the earth’s center. Ascending node is defined as the 
point of intersection of the orbit and the equatorial plane. The orbit passes from south to 
north on the ascending node. Descending node is also defined as the point of intersection 
of the orbit and the equatorial plane. However, the orbit passes from north to south on the 
descending node. Inclination angle is defined as the angle between the orbital plane and 
the equatorial plane. It is denoted by i and measured from the normal of the orbital plane 
to the normal of the equatorial plane. It is the greatest latitude of the orbit, either on the 


southern hemisphere or on the northern hemisphere [12]. 


The argument of perigee is the angle between the ascending node and the perigee. 
It is denoted by y and is measured in the orbital plane at the center of the earth, in the di- 


rection of the satellite motion [12]. 


Right-ascension of the ascending node is measured starting from the line of Aires 
eastward to the ascending node. The longitude of the ascending node is frequently used 
as the right-ascension of the ascending node [12]. The right-ascension is denoted by Q. 
In this study, we assume that the right-ascension is the angle on equatorial plane meas- 


ured eastward from the EO longitude (Greenwich) to the ascending node. 


The mean anomaly is the angular position of the satellite relative to the ascending 


node. The angle is measured from the center of the orbit. This is the angle that a vehicle 
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would travel starting from the ascending node for a given time at a given angular veloc- 


ity. The mean anomaly is denoted by @ [21]. 





\ 
vernal equinox \ 


Figure 14. Orbital elements (Adapted from [20]). 


Kepler’s first law defines the path of the orbit that is followed by a spacecraft as 
an ellipse. The earth sits at the center of the ellipse. The shape and parameters of an ellip- 


tical orbit are shown in Figure 15. 


The semi-major axis is denoted by a and the semi-minor axis is denoted by b [12]. 


The eccentricity is denoted by e and given by Kepler’s first law as [21] 





(3.1) 
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The orbit can be a circle (e = 0), ellipse (0 < e < 1), parabola (e = 1) or hyperbola (e > 1) 
[12]. 





oi we 


Figure 15. The semi-major axis and the semi-minor axis of an ellipse. 


Kepler’s second law assumes that the area on the orbital plane subtended by each 
vehicle’s travel at equal time intervals will be equal to each other [21]. Kepler’s third law 
defines the relationship between the semi-major axis of the orbit a and the angular mean 


motion of the vehicle n, (rad/s) as given by [22] 


gate (3.2) 


where yu, is earth’s geocentric gravitational constant (3.986005e14 m*/s? ). 


The preceding equation assumes a perfect orbit without any perturbation effects. 
This law shows that there is a fixed relationship between period and the size of the orbit: 
the square of the angular mean motion of the spacecraft is inversely proportional to the 
cube of the semi-major axis [12]. 

2. Orbit Perturbations 

The Keplerian orbit assumes that the earth is a perfect sphere with a uniform 
mass, and the gravitational force and the centrifugal force caused by the spacecrafts’ mo- 
tion are the only forces acting on the spacecraft [12]. However, the orbit will suffer some 
disturbance due to several effects. Earth’s oblate shape is the major perturbation that af- 
fects the orbit. There are several other causes of perturbation, such as gravitational forces 
of the sun and the moon, atmospheric drag and solar radiation pressure [22]. The gravita- 


tion forces of the sun and the moon are negligible for the spacecrafts moving on a low 
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orbit. The oblate shape of the earth effects the mean motion of the orbit leading to a new 


definition for the mean motion [12]: 


K, (1-1.5sin’/) 


3.3 
a (1-e ) : re 


n=n,| 1+ 
In (3.3), Ki is a constant (66,063.1704 km). The equation for the nominal mean motion 


i (3.4) 


The perturbation effects other than earth’s oblate shape are neglected in this study. It is 


is derived by re-arranging (3.2) as [12] 





assumed that earth’s oblate shape only affects the mean motion of the vehicle in orbit. 
Cc. ASSUMPTIONS AND REQUIREMENTS 

The hostile ICBMs are assumed to be located in Kilju-kun, North Korea (N 41° E 
129°), Xining, China (N 36° E 101°) and Bushehr, Iran (N 29° E 51°). The ICBMs of 
these rogue states are assumed to have the same characteristics. The mathematical model 
for these ICBMs were introduced in Chapter II. The geographical locations of these bases 
are shown in Figure 16 — 18 [23]. The images of the figures were generated using 
Google Earth® on 5 June 2005. Recently, Iran has been seeking opportunities to improve 
the nuclear reactor in Bushehr. There are some clues to prove that Iran has some missile 


activities in Bushehr [24]. 
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Figure 16. The geographical location of Kilju-kun Missile Base, North Korea (After 
[23]). 





Figure 17. The geographical location of Xining, China (After [23]). 
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Figure 18. The geographical location of Bushehr, Iran (After [23]). 


The ICBMs that are launched from the above sites are assumed to target San 
Francisco, California. The orbit of the EKV carriers should cover all three sites, i.e., the 
inclination angle and the right-ascension angle should meet this requirement. The EK Vs 
that are in the orbit will have the same specifications and will be used against any of the 
three launch sites. Consequently, we should select a circular orbit to guarantee uniform 


coverage (e = 0). 


The intercept should take place after the boost-phase of the ICBM and prior to the 
release of the RVs and decoys. In addition, the EKV is required to kill the ICBM with the 
damage caused by physical impact. In other words, the EKV should achieve a hit-to-kill 
intercept. The time interval (window) for the intercept is assumed to be 15 s, and the in- 
tercept should occur at the end of the boost-phase (three minutes after launch) and before 
the release of the RVs (three minutes and 15 s after launch). In this study, we ensure that 


the intercept occurs before the delivery of the RVs. 


The altitude of the orbit depends on several requirements and limitations. The 
maximum range of the EKV within the intercept window is the major limitation for the 
orbit altitude. The orbit altitude should give enough time and distance to allow the EKV 
to accomplish its mission. That is, if you put the orbit at a very high altitude, the EKV 


will not be able to intercept the ICBM before it delivers the RVs. The capability of the 
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space launcher that will deliver the EKV carriers to the orbit is another important factor 
for determining the altitude of the orbit. The altitude of the orbit is also a measure of the 
lifetime of the system. 
D. ORBIT DERIVATION 

We determine the parameters of a circular orbit according to the assumptions and 
requirements described above. The parameters to be calculated include the inclination 
angle i, altitude of the orbit ,, and the right-ascension of the ascending node Q. 

1. Analytical Solution for Orbit Determination 

Since we selected the circular orbit, the semi-major axis and the semi-minor axis 
are both equal to the sum of the altitude of the orbit and the radius of the earth. The ar- 
gument of perigee, the mean anomaly and the eccentricity are also defined inherently in 
the circular obit (e = 0). After calculating these parameters and knowing the maximum 
range of the EKV, the required number of EKV carriers may be estimated. The parame- 


ters of interest are shown in Figure 19. 
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Figure 19. Orbital plane and the launch sites. 
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a) Inclination Angle 
The inclination angle is determined such that the EKV carrier passes over 
all the launch sites. The orbit of the EKV carriers defines an orbital plane that has the 
earth’s center as the origin and is shown in Figure 19. The equatorial plane includes the 
equator as a circle and has the same origin with the earth. The angle between these two 
planes at the center of the earth is the inclination angle. The ICBM launch sites are de- 
noted by vectors X, (North Korea), X> (China) and X; (Iran) in Figure 19. Ideally, these 
sites would lie on the same global circle to let us use one orbit for all. These vectors 
should be members of the orbital plane by definition. The launch locations are given in 
geodetic coordinates, so we need to convert to the ECEF coordinate system before doing 
further calculations. The locations are first transformed to spherical coordinates and then 
to the ECEF coordinates for convenience. The conversion from spherical coordinates to 
the ECEF is given by [25] 
x R, sin @cos¢ 
X =| y|=| R,sinOsing (3.5) 
Zz R, cosO 


The spherical angles @and ¢ are illustrated in Figure 20. 





Greenwich 


Figure 20. The spherical angles of the given location. 
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The spherical angle @ for the locations in the northern hemisphere is calculated using the 


latitude as 


6-2 66) 


ee 
2 
where wis the geodetic latitude. The spherical angle ¢ is the longitude. The range R, is 
the radius of the earth. 

Kilju-kun, North Korea: 

The spherical coordinate parameters for Kilju-kun are R, = 6371 km, 


0=90° —41° = 49° and ¢=129" . The vector X, can now be calculated in the ECEF coor- 


dinate system using (3.5) as 


-3025.93 km 
X, =| 3736.72 km 
4179.75 km 
Xining, China: 
The spherical coordinate parameters for Xining, China are R,= 6371 km, 
0=90° —36’ = 54’, g=101°.. The vector XY, can now be calculated in the ECEF coordi- 


nate system using (3.5) as 


-983.48 km 
X, =| 5059.55 km 
3744.78 km 


Bushehr, Iran: 


The spherical coordinate parameters for Bushehr, Iran are R= 6371 km, 
0=90° —29° =61’ and ¢=51°. The vector X3 can be calculated in the ECEF coordinate 


system using (3.5) as 


39 


3506.70 km 
X, =| 4330.41 km 
3088.72 km 


The Xi, X2 and_X3 vectors now define the orbital plane and the inclination 
angle. The inclination angle can also be described as the angle between the normal of the 
orbital plane and the normal of the equatorial plane. The normal vector of the orbital 


plane h is calculated by taking the cross product of X, and X; vectors and is given by 


6558382.10 
h=| -24003403.75 (3.7) 
26207085.48 


The normal of the equatorial plane can be any vector along the z-axis, so we can select an 


arbitrary vector on z-axis as 


0 
Z=| 0 (3.8) 
6371 km 


The dot product of the normal of the orbital plane and the normal of the equatorial plane 
divided by the product of the magnitudes of these vectors gives us the cosine of the incli- 


nation angle. Hence the inclination is calculated as 


i=arccos| —_—~—_ 
(atiz | 


) = 0.76 radians or 43.52° (3.9) 
Now we know the inclination angle of the circular orbit. We will target this orbit around 
the earth by calculating the right-ascension angle. 

b) Right-ascension of the Ascending Node 

Right-ascension angle is measured eastward from the longitude 0 (Green- 
wich) to the ascending node of the orbit. In other words, it is the longitude of the ascend- 
ing node. As shown in Figure 19, the line of nodes is the intersection of the equatorial 
plane and the orbital plane. The line of nodes will pass from the ascending node and the 
descending node. The cross product of the normal vectors of these two planes will define 
the line of nodes. The line of nodes ‘PY is given by 
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W=Zxh (3.10) 
which is defined from the center of the earth in the ECEF coordinate system as 


xy 152925685295.96 
Y=! y, |=] 41783452372.92 (3.11) 
Z, 0 


y 


The line of nodes lies along the equatorial plane, which is ensured by the z value of the 
vector being zero. The tangent of the right-ascension angle (.) is defined by the x and y 


elements of the line of nodes vector. The right-ascension angle is computed as given by 


x 


Q= win) = 0.27 radians or 15.28” (3.12) 


Since the right-ascension angle is the longitude of the ascending node, the orbit crosses 
the equator from south to north at longitude 15.28°. 

° Altitude of the Orbit 

The EKV carriers will encounter atmospheric drag in low orbits. This drag 
will increase exponentially as the altitude of the orbit decreases. For example a 100-km 
decrease from 700-km of altitude will increase the drag to be overcome by five times 
[20]. 


The EKV carriers use their onboard thrusters to make the orbital correc- 
tions caused by the atmospheric drag, along with other altitude-dependent effects. After 
the depletion of the thruster fuel, the EK V carriers will start to loose altitude and reenter 
the atmosphere, which will be the end of their lifetime. Since the EKV carriers can carry 
a limited amount of thruster fuel, the altitude of the orbit will determine the lifetime of 


the system [20]. 


The launch vehicle that places the spacecraft into the orbit also has a limi- 
tation for the altitude with a given mass of payload [20]. Since we have a fixed mass of 
EKV, the mass of the EKV carrier is dependent on the booster fuel embedded into the 
EKV. Neglecting the thruster fuel that is needed to keep the EKV carrier in constant or- 
bit, the altitude is a trade off between the mass of the booster fuel on each EKV and the 


capability of the space launch vehicle. As the booster fuel of the EKV increases, its maxi- 
4l 


mum range will increase, which will let us choose a higher orbit. However, the increase 
in the booster fuel will increase the total mass of the EKV carrier, which will decrease the 


orbit altitude due to space launch vehicle payload mass limitations. 


We have selected Titan IV as the representative launch vehicle because it 
is one of the most powerful expendable launch vehicles of the U.S. The total mass of Ti- 
tan IV at launch is 1038 tons [26]. The payload capability of the TitanIV is up to 10,000 


pounds for geosynchronous orbits and is up to 39,100 pounds for low earth orbits [28]. 





Figure 21. Titan IV at liftoff (From [27]). 


The payload system mass as a function of the circular orbit altitude is 
shown in Figure 22. For example, the Titan IV can deliver a 2000 kg of payload into a 


circular orbit with an altitude of 1000 km. 
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Figure 22. Space launcher (Titan IV) payload capacity (After [20]). 


The resulting altitude of the EK V orbit then influences the maximum 
flight range of the EKV in a given time and the altitude of the desired intercept point. We 
assume a 30-s delay for the EKV firing after the [CBM launch, which includes the launch 
detection time, firing decision time and the human interface and communication laten- 
cies. We want to intercept the ICBM before the release of the RVs, which leaves the 
EKV two minutes and 45 s to intercept the ICBM. The maximum range of the EKV as a 


function of the orbit altitude for two minutes and 45 s is shown in Figure 23. 
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Figure 23. Maximum range of the EKV during a maximum allowable period of 2 min- 
utes and 45 seconds 


The maximum range of the EKV decreases as the orbit altitude increases 
because the initial velocity of the EXV decreases as the altitude decreases, which is indi- 
cated by (3.3). 


We used the same motion modeling method described in Chapter II to de- 
termine the maximum range of the EK V. The ICBM launched from Iran has the mini- 
mum altitude at the end of the boost-phase, as shown in Figure 24. This altitude is the 
lower bound, which will require the EKV to travel the maximum range. The altitude of 
this ICBM three minutes and 15 s after the launch is 213 km. The geometry of the trajec- 
tories of all three ICBMs are shown in Figure 7 in Chapter II. 
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Figure 24. The altitudes of the ICBMs at the end of boost-phase. 


The altitude of the orbit and the maximum range of the EKV together de- 
termine the down-look launch angle & and the coverage angle @ as shown in Figure 25. 
Here we assume that the EKV cannot achieve a successful intercept after it passes over 
the desired intercept point due to its huge initial velocity in the direction of the orbit. 
Maximum Probable 
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Figure 25. Orbital plane and the intercept geometry. 
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The maximum probable orbit altitude is determined by the mass of the 
payload and the capacity of Titan IV together, as shown in Figure 22. The down-look 
launch angle & is calculated using the geometry in Figure 25 and is shown in Figure 26 as 


a function of orbit altitude. 
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Figure 26. The down-look launch angle as a function of the orbit altitude. 


The coverage angle @ is measured at the center of the earth from the EKV posi- 
tion vector to the desired intercept point as shown in Figure 25. The coverage angle is 
calculated using the geometry in Figure 25 and is shown in Figure 27 as a function of or- 


bit altitude. 
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Figure 27. Maximum coverage angle as a function of orbit altitude. 


The coverage angle of the EKV decreases as the altitude increases, which 
makes us choose a lower altitude orbit. The total mass of the system is 2000 kg, which 
dictates that the maximum orbit altitude be 1000 km due to space launch vehicle limita- 
tions, as shown in Figure 22. For an orbit altitude of 1000 km, a coverage angle of ap- 
proximately 10.6 degrees results, as shown in Figure 27. Consequently, to cover 360 de- 


grees requires [360/10.6] =34 launchers. 


Now we defined the required orbit and its elements as a circular orbit with 
an altitude of r,= 1000 km, an inclination angle of i= 43.52° and a right-ascension angle 
of Q = 15.28". Since the orbit is circular, the eccentricity is approximately zero and the 
semi-major axis and semi-minor axis are both equal to the sum of the radius of the earth 
and the altitude of the orbit. Next we will derive the orbital parameters in the ECEF coor- 
dinate system and implement this orbit in the model we developed. 

2. Numerical Implementation of the Selected Orbit 
Since the orbit is defined in its own specific elements, we need to define it in the 
ECEF coordinate system to integrate with the model we have developed. First we will 


calculate spherical coordinates of the orbit and then transform them into the ECEF coor- 
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dinate system. The spherical coordinate angles, the ECEF coordinate system and the orbit 


parameters are illustrated in Figure 28. 
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Figure 28. The orbit parameters in the ECEF and spherical coordinate systems. 


In Figure 28, Xzxr is the current position of the EKV, is the right-ascension 


of the ascending node and y, is measured from the projection of Xzxv on the equatorial 


plane to the line of nodes. The length of the Yzxr is the sum of the orbit altitude and the 


radius of the earth. The spherical angles @and ¢are also shown in Figure 28. The mean 
anomaly @ is measured from the ascending node at a given time. The time starts when the 


EKV crosses the ascending node and the mean anomaly as a function of time is given by 
@=nt (3.13) 


After some geometry, we calculated the spherical angles as a function of inclina- 


tion angle and the mean anomaly and tabulated them in Table 6. 
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Table 6. Spherical coordinates calculation for EKV at a given time. 
Spherical Coordinates 
@ (rad) @ (rad) r=7,tK, 
0<<F Fo eos! (Veos® «+ sin’ «o-c0s") Q+ tan“ (tan @-cosi) 7371 km 
F<ose F008" (eos? + sin® cos") Q+2-tan'(-tan@-cosi) | 7371 km 
n<os3e F+e0s"( cos’ @+sin* @-008"/) Q+n+tan'(tanw-cosi) | 7371 km 
at< o<2n F +e0s"'( cos’ @+sin" @-008"/) Q4+2x-tan'(-tan@-cosi) | 7371 km 

















After calculating the spherical coordinates, the EKV position vector is converted 


to the ECEF coordinate system as given by [25] 


The above calculations are implemented and simulated in MATLAB? with the given 


rsin@cos@ 


X xy =| rsin@sing 


rcos@ 


ECEF 


(3.14) 


conditions of @. The orbit and the target ICBM locations are shown in Figure 29. This 


numerical computation of the orbit is used later to calculate the initial launch parameters 


of the EKV against the ICBMs that are launched from North Korea, China or Iran. 
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Figure 29. The orbit of the EKV and the target ICBM locations. 


E. SUMMARY 

In this chapter, an orbit is derived to cover the three example ICBM launch sites 
uniformly and continuously. The lifetime of the space system, the capabilities of the EKV 
and the space launch vehicle (Titan IV) are taken into account in deriving the required 


orbit. After selecting the orbit, it is simulated in MATLAB". 


50 


IV. EKV GUIDANCE METHODS 


This chapter will investigate different guidance rules for space-based interception 
of hostile ICBMs launched from North Korea, China or Iran. The pursuit guidance (PG), 
proportional navigation guidance (PNG) and bang-bang guidance (BBG) rules will be in- 
troduced and a three-dimensional implementation will be presented. The total com- 
manded acceleration, intercept time, and miss distance generated by these guidance rules 
for the given scenario will be the major parameters considered to select one of them or 
combine them properly. 

A. DESCRIPTION OF THE SCENARIO 

The interceptor missile is a one-stage, boosted, exo-atmospheric kill vehicle 
(EKYV) orbiting in a circular LEO orbit with an inclination angle of i = 43.52’, a right- 
ascension angle of Q =15.28° and an altitude of r.=1000 km. The target ICBMs are com- 
ing from North Korea, China and Iran. The coordinates of these launch sites are N49- 
E129, N36-E101 and N29-E51, respectively. The ICBMs are targeting the city of San 
Francisco, California. 

B. COORDINATE SYSTEMS AND COORDINATE CONVERSIONS 

In this chapter, three different coordinate systems will be introduced. First, the 
common coordinate system that will be used in the calculations will be the Earth-centered 
Earth-fixed (ECEF) coordinate system (Xg,Ye,Zr). The ECEF coordinate system is a 
three-dimensional orthogonal Cartesian coordinate system with the origin at earth’s cen- 
ter. The x-axis passes through Greenwich (E0), the y-axis passes through E90 and the z- 
axis passes through the North Pole. The spherical earth rotates about z-axis. We assume 
that the earth is a perfect sphere and the angular rotation of the earth will be involved in 
the computations. The second coordinate system is the Airframe Body Coordinates 
(ABC), which is a rotating system with the airframe of the vehicle (Xm, Ym, Zm for EKV 
or Xr, Yr, Zr for ICBM). The x-axis of this coordinate system lies within the velocity 
vector of the vehicle, assuming that the velocity vector and the missile body are aligned. 
The y-axis is the left wing of this vehicle and the z-axis is orthogonal to the x and y-axes 
complying with the right hand rule [29]. The third coordinate system is the line of sight 
coordinate system (LOS). This is a two-dimensional coordinate system, which can also 
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be referred to as the LOS plane (Y,, Z,). The EKV position relative to the ICBM is de- 
fined by using this coordinate system. The distance Y, is the LOS vector projected on the 
equatorial plane and the distance Z, is the LOS vector projected on the z-axis of the 
ECEF coordinate system [30]. The LOS plane and the ECEF coordinate system are illus- 


trated in Figure 30. 
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Figure 30. Guidance coordinates geometry: ECEF coordinate system, LOS plane, ABC 
system. 
The missile and target parameters, such as velocity and acceleration, are defined 
in the ABC coordinate system, and the computations are conducted in the common coor- 
dinate system (ECEF). Hence, we need to define and calculate the coordinate conversions 


to transform all the coordinates to the ECEF coordinate system. 
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The acceleration command will be applied to the yaw and pitch axes of the EKV, 
which are represented by Xm and Yw, respectively. In order to represent these accelera- 
tion forces in the ECEF coordinate system, we have to rotate the coordinate system about 
the rotating axis. The angles resulting from these rotations are called the Euler angles 


[31]. The rotation matrix about the x- axis by angle y can be written as 


1 0 0 
Y=/0 cosy sinw (4.1) 
|0 —siny cosy 


The rotation matrix about the y-axis by the angle @can be written as 


[cos@ 0 —sin@ 
O=| 0 1 0 (4.2) 
[sind 0 cosé 


The rotation matrix about the z-axis by the angle ¢ can be written as 


[ cosd sing 0 
@M=|-sing cos¢? 0 (4.3) 
0 0 1 


For consecutive rotations, we need to multiply the corresponding rotation matrices [32]. 
In order to transform the ECEF coordinate system into ABC or LOS coordinate systems, 


we first rotate the EK V about the y-axis by the pitch angle 0. The next necessary rotation 





is about the z-axis by yaw angle g. The required overall rotation matrix can be generated 
by multiplying these two rotation matrices as 

C=00 (4.4) 
After some matrix algebra and using appropriate angles, we can write the open form of 


general conversion matrix C as follows 


cosOcosé cos@sing sind 
C= —sing cos¢ 0 (4.5) 
—sin@cos¢ —sin@sing cosé 
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Using (4.5), we can now derive the transformation vectors for the defined coordinate sys- 
tems. First, we define coordinate transformation from the ECEF coordinate system to the 
LOS coordinate system as [30] 

i, i, 

J j= G| a (4.6) 

k, k 


L 1 


where i, j and k represent the unit vectors in x, y and z axes, respectively, in the corre- 
sponding coordinate system. The conversion matrix C, is defined as a function of 

6, and ¢, angles, which are shown in Figure 30 and represent the appropriate Euler an- 
gles [30]. The subscripts /, L and M represent the ECEF coordinate system, the LOS co- 
ordinate system and the ABC system, respectively. The coordinate conversion from the 


LOS coordinate system to the ABC system using (4.5) is defined as 


Ju J=Cn| Jt (4.7) 


The conversion matrix C,, is defined as a function of @, and ¢, angles, which are shown 
in Figure 30. The 6,, and¢,, angles are the spherical coordinate angles for the velocity 


vector of the missile: 


8, =O, -9, (4.8) 
$, = Pu 4% : 


Cc. PURSUIT GUIDANCE 

The aim of the pursuit guidance, also referred to as beam rider guidance, is to 
keep the EKV velocity vector on the target and to finally hit the target. The pursuit guid- 
ance is very simple to implement and does not require a complicated seeker design. It re- 
quires a velocity and geometry advantage to perform a successful intercept. This guid- 
ance rule performs better in tail aspect intercepts with a velocity overtake. It has some de- 


ficiencies for cross-range collision geometries [14]. 
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1. Principals of Pursuit Guidance 

The seeker of the missile is required to produce only a line of sight angle for the 
pursuit guidance. Since the EKV has very efficient infrared sensors on board, the line of 
sight angle can be measured very accurately in the end game. Pursuit guidance also re- 
quires the measurement of the angle of the missile velocity vector. Both the line of sight 
angle and the missile velocity vector angle should be referred to the same fixed plane, 
which is the equatorial plane in this study. The seeker measures the line of sight angle in 
the ABC and is then converted to the ECEF coordinate system. The inertia navigation 
system (INS) will also measure the missile velocity vector angle in ABC and is then con- 


verted to the ECEF coordinate system. 


The angle error is generated by taking the difference between the LOS angle and 
missile heading, as shown in Figure 31. After capturing the angle error, the acceleration 
command is generated from this error and guidance gain K. The command acceleration is 
filtered by the 7(s) and the missile dynamics uses this acceleration to change the course 
of the missile by applying the thrusters. Since the interception will take place in space, 
the guidance maneuver will be accomplished using thrusters instead of aerodynamic sur- 


faces. 
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Figure 31. Pursuit guidance flow chart. 


2. Three-Dimensional Implementation of the Pursuit Guidance 

Since the IR seeker of the EKV measures angles very accurately, we can assume 
an exact measurement of the target position. For the pursuit guidance, the seeker is re- 
quired to provide the guidance system with the LOS angle. The seeker will measure the 
yaw and pitch angles of the EKV relative to the ICBM within its spherical ABC system. 


These yaw and pitch angles are represented by ¢,, and 0, 


respectively (see (4.8)). Be- 


fore defining the yaw and pitch angles, we first need to define the LOS vector. The LOS 
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vector can be calculated by vector subtraction of the target and missile position vectors. 


The target and missile state vectors are defined as 


Sel yn 4% 4) 


; ; 7s (4.9) 
x u> [Xn Yn 2m Xm Vn Zn] 
Using (4.9), we can define the LOS vector as 
% Xm 
X= X;,(1:3)-Xy (1:3) =] ¥, |=] In (4.10) 
2; 2, ma 2m 
The line of sight (LOS) angles 6, and ¢, are derived as given by 
(4.11) 





The missile heading is computed by using the velocity vector of the missile, 
which is defined in the state vector, assuming that the missile body coincides with the ve- 
locity vector. The velocity vector contains the last three components of the state vector. 


The missile heading angles @,, and ¢,, are computed by 


(4.12) 





Now we can define the yaw and pitch angles as a function of LOS angles and the missile 


heading. The yaw and pitch angles are shown in Figure 30. 


These 6,, and ¢,, angles represent the missile heading error. Pursuit guidance at- 
tempts to zero this error that will guide the EKV to the ICBM. The acceleration command 
is produced for both yaw and pitch axes to correct this angle error. The yaw and pitch an- 


gle deflections are converted to ECEF components as [33] 
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A, =, Sin 0, 
oo) (4.13) 


4, = $, cos 4, 


where / represents the angular error. The guidance commands are applied to be perpen- 
dicular to the LOS vector. However, we can only apply acceleration commands perpen- 
dicular to the EKV body because of fixed thruster orientation. Hence, using the equations 
(4.6) thru (4.7), we define the yaw and pitch turn rate commands for pursuit guidance as 


(After [33]) 


Oyu = KA, sin 8, sing, —KA, cos8,, (rad/s) 


(4.14) 
Oynen = KA, COS g,, (rad/s) 


where K is the guidance gain. Increasing the guidance gain will increase the turn rate, 
which will make the guidance more responsive. However, it will also make the guidance 
more susceptible to saturation. Hence, the selection of the guidance gain is a trade off be- 


tween agility and stability. 


Definition (4.14) yields the turn rate commands, but we still need the acceleration 
command to apply to the thrusters of the EK V. The acceleration command is evaluated 


using coordinated turn definitions, which is sketched in Figure 32 [19]. 


v 


Figure 32. Coordinated turn schematic (After [34]). 
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The acceleration is defined as a function of velocity and the radius of the turn circle [15]: 
Vy? 
a=— (m/s’) (4.15) 
r 
where a is the acceleration that is normal to the velocity vector, V is the tangential veloc- 
ity and r is the radius of the turn circle. The turn rate can be defined as the total angle 
traveled over the time. The total angle traveled in a complete turn is 27 radians and the 


time to travel this angle is 





2ar 
t=—— (s 4.16 
= (s) (4.16) 
Hence, we can write the turn rate as 
2a V 
o= =— (rad/s 4.17 
2ar/Veor ( i ) ( ) 


Substituting (4.17) into (4.15), we can define the relationship between the acceleration 


and the turn rate as 
a= (m/s’) (4.18) 
Now we can write acceleration commands instead of turn rate commands as (After [30]) 


Nyy, = KVyA, sin 8, sing, —KV,,2, cos, (m/s”) 


ce (4.19) 


N joy = KV,,A, COS$,, (m/s’) 


pitch 


The yaw and pitch acceleration commands should be expressed in the ECEF co- 
ordinate system to be implemented in the simulation. Using the unit vector representa- 


tion, we write the resultant acceleration command in ABC as 


tl 


Ayan =[9 Mya Ryne |] Ju (4.20) 
ky 
where n is the acceleration command vector in ABC. Using (4.7), we convert this 


pursuit 


acceleration command into LOS coordinates as follows 
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i, 
Mpursit =[9 yaw Mncr |On| Je (4.21) 
k, 


Finally, substituting (4.6) into (4.21), we define the acceleration command in the ECEF 
coordinate system as 
i, 
Myatt =[9 Brow pace |nEx | Jr (4.22) 
k, 


Latencies in the seeker measurements and latencies in applying the acceleration 
command are two dominant factors that are affecting the intercept capability of the guid- 
ance. These two factors are modeled in this study by using a third order transfer function. 
These latencies will cause a significant miss distance, which would not take place if there 


were no latencies [35]. The transfer function is defined in the Laplace domain as [9] 


7(s)-__ | (4.23) 


(5) 


where n is the order of the transfer function and 7 is the time constant in seconds. The 
time constant can vary from 0.5 s to 2 s depending on the agility of the missile [15]. The 
intercept geometry for space-based intercept of an ICBM generally has a head-on aspect. 
This geometry results in very high closing velocities and requires a highly responsive 
EKV. The EKV that is modeled in this study is a hypervelocity advanced missile with a 


time constant of 0.5 s. The transfer function for c = 0.5 scan be derived as 


216 


T9)=3 ODS 
) 5° +1857 +1085 +216 


(4.24) 


The time domain response is derived by taking the inverse Laplace transform of 7(s) and 


is given by 





(= 2"{ 7) 1-80 sorrow (4.25) 
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The time domain response is computed using MATLAB” and the resulting plot is 


shown in Figure 33. 
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Figure 33. Time domain response of the EKV guidance system. 


As arule of thumb, four times the time constant shows us when the transition will 
end and the command will be fully applied [15]. Here the time constant is 0.5 s, so the 
transition should die out in 2 s. As shown in Figure 33, the time domain response reaches 
0.9995 at 2 s, which can be considered the end of the transition. 

D. PROPORTIONAL NAVIGATION GUIDANCE 
The Proportional Navigation Guidance (PNG) produces a perpendicular accelera- 


tion that is a function of the LOS angle rate 6, , the closing velocity V. and the propor- 


tional navigation constant N. The turn rate command is generated proportional to the an- 
gular velocity of the LOS [36]. The EKV seeker provides an accurate LOS angular veloc- 
ity measurement, which will increase the performance of the intercept in our model. The 
idea is to guide the EKV directly to the collision point by leading the target, so the EKV 
can intercept the ICBM at an optimum time and with minimum command effort. The 
PNG is well known for minimizing the total control effort. Theoretically, the guidance 
command should be applied perpendicular to the LOS. However, in practice we can only 


apply the guidance command perpendicular to the missile body [33]. The guidance com- 
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mand derived below assumes that EK V is launched towards the missile, but we will see 


that it is still optimal when the scenario does not meet this requirement. 


First, we will discuss the PNG in two-dimensions to have a better understanding 
of the theory. The scenario is defined in Figure 34 for the two-dimensional geometry. The 
geometry assumes a constant speed impulsive target and a constant speed impulsive in- 
terceptor. The collision flight path is defined from the interceptor’s initial launch point to 
the intercept point. Note that the LOS angle remains unchanged on this collision flight 
path. The PNG aims to keep the missile on this path by producing acceleration com- 
mands. These acceleration commands should also assure that the LOS angular velocity 
stays zero. For the scenario given above, the PNG is the optimal guidance. It will become 
sub-optimal with an accelerating target or an accelerating interceptor. PNG can still be ef- 


fective enough by varying the guidance coefficient, which is discussed below. 


“ 
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Figure 34. PNG collision geometry. 


The LOS angular rate 6, is measured by the seeker. The closing velocity V. is the 
inverse of the time derivative of the distance between the interceptor and the target R. 
The PNG acceleration command perpendicular to LOS is [9] 

n, = NV.O, (4.26) 
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Note that the acceleration command derived above is perpendicular to the LOS 
vector. However, we can only apply command forces perpendicular to the missile body. 
The component of this command perpendicular to the missile body is derived assuming 


that the missile velocity vector is aligned with the missile body [37] as follows 


n NYG. (4.27) 
0s, 
1. Closing Velocity Approach 


The coordinate conversion definitions (4.1) through (4.8) also hold for the PNG, 


and the geometry is defined in Figure 30. The LOS angular velocity vector ©, in three- 


dimensional geometry is defined by using the coordinate conversion vectors [30] 
a=[A A Ad A (4.28) 


where / is the LOS angular velocity in the corresponding component of the LOS coor- 
dinate system [30]. We defined the three components of the LOS angular velocity as 

i, = 4, sin 6, 
i j (4.29) 





A, = 4, 0086, 
The Euler angles used in (4.29) are presented in Figure 30, and the LOS angular veloci- 


ties of 6, and 4, are in the spherical ECEF coordinate system. The seeker of the EKV is 


capable of measuring these angles and passing them to the guidance system. The spheri- 
cal ECEF angles of the LOS vector are derived using (4.11). The angle rates are then 


given as 


(4.30) 





where x,, y, and z, are defined in (4.10). 


We derive the guidance command perpendicular to the EKV as a function of the 


closing velocity, angular rates and the proportional navigation constant N as [9], [30] 


Npxenv =-NV,A, sin 8, sing, + NV,A, cos 8, 


(431) 
NpnGpitch = —NV.A, cos@,, 


where N is the guidance coefficient varying from three to five, depending on the maneu- 
verability of the interceptor and the target. We should select a higher guidance coefficient 
for more agile targets and missiles [9]. The closing velocity is denoted by V,. in the pre- 
ceding equation and is defined as the time derivative of the distance between the intercep- 
tor and the target: 


Vv. =-k= 2 fey +2} (4.32) 


t 


The acceleration command should be derived in ECEF Cartesian coordinates in 
order to implement this theory in the simulation. The acceleration command is first calcu- 
lated in ABC as given by 

iy 
Mev = [0 Merenm  Mencpner I) J (4.33) 
k, 


M 


where npyg is the acceleration command vector in ABC. Using the method (4.20) 


through (4.22), the PNG acceleration command is transformed to the ECEF coordinate 
system: 
i, 
Neng =[0 "pnGyaw Tibicsoey.| Caer di (4.34) 
k, 


The block diagram of the PNG implementation is shown in Figure 35. 
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Figure 35. Block diagram of the PNG. 


2. Missile Velocity Approach 

The angular velocity of the LOS vector is derived using equations (4.28) through 
(4.30) for missile velocity approach. Blakelock [36] states that the acceleration command 
generated by the PNG should be proportional to the missile velocity instead of the closing 
velocity. Blakelock also states hat the PNG gain should be greater or equal to one to en- 
sure an intercept. He recommends selecting a PNG constant between two and six for ac- 
celerating or maneuvering targets [36]. We will denote the PNG constant with N' to 


make a distinction between the closing velocity approach. 


The acceleration commands for the yaw and the pitch axes of the interceptor is 


defined in missile body coordinates as [33] 


Npxouw =—N'V,,A, sin 8, sin g,, + N'V,,A. cos 6, 


aie " (4.35) 
Npyepien = —N'V,,,A, COSY, 


where V, is the missile velocity, which is measured by the onboard INS or the GPS in the 


EKV. The acceleration command vector in the ABC system is obtained as [33] 


ing 
Nyy =[0 May Mpner |] Jar (4.36) 
ky 


The acceleration command is transformed into the ECEF coordinate system using the 
method introduced by (4.20) through (4.22): 


t 


Mpyg =[0 Mow Mpuch [CnC di (4.37) 
k, 
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The block diagram for the PNG law for missile velocity approach is also repre- 
sented by Figure 35, with N replaced by N’ and V.. replaced by Vin. 
E. BANG-BANG GUIDANCE 

Bang-bang guidance (BBG) is a derivative of proportional guidance. Bang-bang 
guidance applies the maximum possible acceleration in the direction of the LOS angular 
velocity. The missiles with thrusters as the control elements apply this guidance very ef- 
fectively. The guidance command is defined as a function of the LOS angle, the LOS an- 


gle rate and the closing velocity [37] 
(4.38) 


where a,, is the maximum applicable lateral acceleration of the missile. 


The BBG produces the maximum maneuver resulting in a large amount of drag. 
This drag degrades the performance of the interceptor in the atmosphere. Since the EKV 


is moving in high atmosphere, this drag will not effect the performance of the EKV [38]. 


Defining the bang-bang acceleration command using the existing derivation for 
the PNG will make this guidance easier to implement into the model that we developed. 
The direction of the guidance command perpendicular to the missile body is derived with 
the same method we used in PNG. The unit vector, which defines the direction of the 
bang-bang acceleration command, is the same unit vector of the PNG guidance command 


and is given by 


fy = oe (4.39) 


Pex 








where fi, is the unit vector of the bang-bang acceleration command. The bang-bang ac- 


celeration command is derived using the preceding unit vector as 
Nyy = AyMyy (4.40) 


Equation (4.39) and (4.40) together define the bang-bang acceleration command 


along the lines of the same derivation for PNG. 
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F. SUMMARY 

In this chapter, pursuit guidance, proportional navigation guidance and the bang- 
bang guidance rules are introduced, and a three-dimensional implementation is mathe- 
matically modeled. We will compare their performance in the next chapter. Also, a hy- 
brid guidance algorithm that combines the principles of these rules is identified in the 


next chapter. 
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Vv. COMPARISON OF GUIDANCE LAWS 


In this chapter the three-dimensional SIMULINK® model is developed and the 
simulation is tested for different scenarios. The performance of the different guidance 
tules are tested for the North Korea case and then an hybrid guidance (HG) algorithm is 
developed based on these results. Here the launch delay is assumed to be zero. The hy- 
brid guidance algorithm is then tested against three example ICBM threats with four dif- 
ferent EKV launch points for each threat. For all guidance simulations, the ideal seeker is 
assumed to see the ICBM at all times. 

A. SIMULINK® SIMULATION MODEL DESCRIPTION 

The ICBM dynamics modeled in Chapter II, the EKV dynamics and orbit parame- 
ters modeled in Chapter III and the three guidance rules developed in Chapter IV are im- 
plemented ina SIMULINK® model. 

1. Simulation Initialization 

The SIMULINK® model requires initial parameters of the EKV and the ICBM to 
run the simulation. These parameters are entered by the user by running SimulationInit.m 
MATLAB* file separately. The initial parameters include the EKV launch point on the 
orbit, the ICBM launch site and the launch delay. 


The initial state vector is generated in the code for the initial point of the EKV on 
the orbit at launch. The user is asked to enter the mean anomaly time, which starts when 
the EKV passes by the equator from south to north, i.e., time zero for the EKV position 
starts at the equator. The entered time is used to calculate the position of the EK V on the 
orbit. The state vector of the EKV is then calculated and passed to the SIMULINK® 


model. 


The launch location of the ICBM and the EKV launch delay are also selected by 
the user. The three choices of launch locations are North Korea, China or Iran. The 
launch location and the initial launch parameters for a San Francisco attack of the se- 
lected ICBM are predetermined. With the parameters above, the initial state vector of the 


ICBM is calculated and passed to the SIMULINK® model. The propellant masses of the 
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selected ICBM are also predetermined. The MATLAB” functions that are used in this 


model are listed in Table 7. 


Table 7. Code listing for the SIMULINK® model. 





















































File Name SIMULINK® Block Comments 
SimulationInitm | N/A Should be run before simulation 
ICBMmotion.m ICBM Dynamics Calculates change in ICBM state 
Seeker.m Seeker Tracks target for LOS rate and V. 
HitDet.m Seeker Detects hit or miss and stops simulation 
PG.m Guidance Calculates PG acceleration command 
PNG.m Guidance Calculates PNG acceleration command 
HG.m Guidance Calculates HG acceleration command 
Limiter.m Guidance Limits the acceleration command 
Mag.m Guidance Calculates the magnitude of a given vector 
EKVmotion.m EKYV Dynamics Calculates change in EKV state 
PlotSimResults.m | N/A Plots the results that are stored in workspace 

2. SIMULINK® Model Description 


The model incorporates ICBM dynamics, seeker, guidance unit and EKV dynam- 
ics in four different subsystems. The overall model is shown in Figure 36. The flow dia- 


gram for the model is provided in Appendix A. 
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WARNING 
This model requires initial parameters 
Initial parameters must be generated by running Simulationinitm seperately 
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or 
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Recommended Configuration Parameters 
(Solver Options) 


Max step size: 0.001 


Min step size: 0.0007 
Solver: ode45(Dormand-Prince) 
Relstive tolerance:1e-6 





Figure 36. Overall space-based intercept SIMULINK® model for PG, PNG, BBG and 
HG against North Korean, Chinese and Iranian ICBMs. 


This model is capable of running for PG, PNG or HG one at a time. Adding new 
guidance methods is also possible with a new MATLAB® script file for new guidance 
methods. The model gives three options as the ICBM launch point. The subsystems of the 
model are ICBM Dynamics, Seeker, Guidance Unit and the EKV Dynamics. 

3. ICBM Dynamics Subsystem 

The ICBM dynamics subsystem reads the required initial parameters from the 
lookup tables. The tables read these parameters from the workspace resulting from Simu- 
JationInit.m file. Initial parameters include the selected launch site, the ICBM initial pro- 


pellant mass and the ICBM initial state vector. 


The subsystem for ICBM dynamics is shown in Figure 37. The JCBMmotion.m 


file calculates the time derivative of the state vector and the propellant mass, and the in- 
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tegrators within the model integrates these vectors and returns the results back for the 


next iteration. 
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Figure 37. SIMULINK™ model for ICBM Dynamics. 


4, Seeker Subsystem 

The seeker provides the guidance unit with ICBM parameters, such as LOS angle, 
LOS angle rate, off bore sight angles, range and closing velocity. The miss distance is 
also measured in this unit, and the simulation is terminated upon a hit or miss. The model 


is shown in Figure 38. 
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Figure 38. SIMULINK™ model for Seeker design of the EKV. 
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5. Guidance Subsystem 

The guidance subsystem of the model takes the seeker outputs as the input and 
uses them to generate the guidance acceleration command. The generated guidance 
command is filtered by a limiter for maximum acceleration capability of the EKV, and 
the total system delay is applied by the autopilot 7(s). The guidance subsystem model is 
shown in Figure 39. The EKV parameters of interest are also stored in the workspace 


variables in this subsystem. 
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Figure 39. SIMULINK™ model for Guidance Unit of the EKV. 


6. EKV Dynamics Subsystem 

EKV Dynamics uses the same method as in ICBM Dynamics. The initial parame- 
ters are read from the workspace by the lookup tables. These initial parameters include 
the initial propellant mass and the initial EXV state vector. The subsystem takes the guid- 
ance command input from the guidance unit and applies it to EKV Dynamics. The EKV 
dynamics subsystem moves and steers the EK V towards the ICBM during its flight and 
returns an EKV state vector as the output. This output is also carried to the seeker subsys- 


tem by a feedback loop in order to model the INS/GPS unit of the EKV. 
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Figure 40. SIMULINK® model for EKV Dynamics. 


The flow chart of the SIMULINK® model and the MATLAB” functions are pro- 
vided in Appendix A. 
B. SIMULATION RESULTS FOR PURSUIT GUIDANCE 

For all guidance laws, the intercept is tested for North Korea in order to develop 
an optimum algorithm. This algorithm is then tested for all launch sites of interest. The 
SIMULINK™ model is set for pursuit guidance against a North Korean ICBM. The accel- 
eration limiter is set to 500 m/s” . The altitude of the EKV is 1000 km and has a maxi- 
mum range of 1520 km (see Figure 23). With these parameters, it is predicted that if the 
EKYV lies within a coverage angle of 10.6° about the launch site (see Figure 25 and 
Figure 27) it will be able to hit the ICBM successfully. This maximum coverage con- 
straint will be met if the EKV is launched where the mean anomaly (angle from equator) 
of the orbit is 98.6’. This angle is the position of the EKV carrier on the selected orbit. 
This position of the EKV carrier corresponds to 1725 s after its pass by the ascending 
node and lies at the edge of the coverage region. The initial velocity of the EKV is deter- 


mined by the angular velocity of the EKV carrier on the orbit (7.35 km/s) and it is tan- 
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gent to the orbit at the given launch time. The initial (t = 0) state vector of the EKV is 


calculated for the parameters above as 
Sig -2456646.97 


yy | | 4807702.67 
zy |_| 5018218.97 


X,,(0)= = 4.41 
w (0) ky -6436.30 a6)) 
dy 876.06 
iy 3447.19 | cre 


The initial state vector of the target ICBM is defined by its launch location and 
the silo exit velocity. The silo exit velocity of the ICBM is 18 m/s[6]. The location of the 
North Korean ICBM is N41-E129 (see Table 3 for initial launch parameters). Using the 


required conversions, the initial (t = 0) state vector of the ICBM is obtained as 


Xp -3025932.75 
Vr 3736715.75 


z, |_| 4179752.07 
X,(0)=| 7 |= (4.42) 
ky -8.76 
vp 9.48 
a 12.54 |e 


The EKV and ICBM dynamics are simulated as discussed in Chapter II. The 
ICBM is a three-stage missile targeting San Francisco, California. The ICBM reaches a 
velocity of ~6.5 km/sand an altitude of 182 km of altitude at the end of the boost-phase 
(three minutes). The EKV is a single stage kill vehicle with control thrusters that has an 
initial velocity of ~7.35 km/s and altitude of 1000 km. The EKV guidance coefficient K 
is set to 3, which will introduce the least amount of oscillations in the acceleration com- 
mand while keeping the EKV agile. 

The three-dimensional encounter geometry is shown in Figure 41. It is shown in 
the geometry that the EKV falls behind the ICBM, which is the nature of the pursuit 
guidance. After falling behind the ICBM, the EKV uses its velocity advantage to over- 
take and catch the ICBM at the intercept point. 
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Space Launch Point 





Figure 41. The intercept geometry for PG with K = 3. 


Figure 42 shows the miss distance for the PG intercept scenario as a function of 
time with K = 3. The minimum miss distance of this interception is 0.04 m. The mini- 
mum miss distance is small enough to consider this intercept as a hard-kill-hit because 
the diameter of the ICBM is 2.3 m [6]. Figure 44 shows the altitude of the EKV as a 
function of scenario time for PG with K = 3. The EKV falls behind the ICBM during the 
scenario and can be observed in Figure 43, where the EKV altitude is below the ICBM 
altitude. The interception time is read from Figure 42 as 3.93 minutes, which is greater 
than our assumption, which was 3.0 to 3.25 minutes, i.e., the RVs have already been 
launched prior to the intercept. This latency is mostly caused by the lagging trajectory of 
the EKV, which causes it to travel a greater distance than expected. The altitude of the 
EKV increases at the beginning of the scenario because it does not point to the ICBM at 
the launch time, but the pursuit guidance succeeds in providing a positive closing velocity 


after ~0.3 minutes, and the EKV starts descending toward the ICBM. 
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Figure 42. Miss distance for pursuit guid- Figure 43. Altitude of the EKV and the 
ance with K = 3. ICBM for pursuit guidance with K = 
3. 


The performance of the pursuit guidance is not satisfactory mainly due to the ex- 
tra distance the EKV travels under PG. Figure 44 shows the instantaneous acceleration 
applied by the EKV under PG with K = 3. Figure 45 shows the cumulative acceleration. 


The instantaneous acceleration command reaches the maximum applicable value 

(500 m/s’) at the beginning of the interception in order to point the EK V towards the 
ICBM, as shown in Figure 44. This maximum acceleration continues for 30 s until the 
EKV starts pointing towards the ICBM. The total acceleration command effort applied to 
the EKV can be read from Figure 45 as 33,000 m/s°. This total acceleration command 
gives an idea about the total effort used by the EKV to hit the ICBM. Minimizing this to- 
tal command effort will decrease the cost proportionally [39]. The EKV begins to fall be- 
hind the ICBM ~130 s after the launch and starts to apply more acceleration command to 
make the EKV point to the ICBM again. Here, acceleration command increases gradually 
and reaches 212 m/ s° . The total acceleration command increases linearly up to 30 s, 
when the maximum acceleration command is applied continuously, as shown in Figure 
45. The oscillations observed in the instantaneous acceleration are the result of the auto- 


pilot 7(s). 
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Figure 44. Instantaneous acceleration ap- Figure 45. Cumulative acceleration applied 
plied on the EKV for pursuit guid- on the EKV for pursuit guidance 
ance with K = 3. with K = 3. 


Figure 46 shows the speed of the EK V under PG and the speed of the ICBM. 
Figure 47 shows the closing velocity. The initial velocity of the EKV is a result of the or- 
bital angular velocity and gives an advantage to the EKV for this scenario, where the ini- 
tial velocity points in the direction of the ICBM within ~60°. The velocity of the EKV 
starts to decrease, as shown in Figure 46, due to gravitational forces. After the EKV starts 
to point to the ICBM, the velocity of the EKV starts to increase again due to the booster 
thrust. The rapid increase in velocity stops when the booster fuel is expended, but the in- 
crease in velocity is sustained due to gravitational forces. The EKV falls behind the 
ICBM ~130 s after the launch because of the nature of the pursuit guidance. Hence, the 
velocity of the EKV starts to decay because it starts to move against the gravitational 


forces after this point, while the ICBM continues accelerating. 


The closing velocity of the EKV to the ICBM is negative at the beginning because 
the initial velocity of the EKV is not pointing to the ICBM due to the orbital angular ve- 
locity vector at launch. The pursuit guidance commands the missile to turn towards the 
ICBM, which results in a rapid increase in the closing velocity (as shown in Figure 47) 
because of the favorable interception geometry. The closing velocity gradually increases 
until the EKV falls behind the ICBM and rapidly decreases until the end of the boost- 
phase of the ICBM. 
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Figure 46. Speeds of the EKV under PG Figure 47. Closing velocity of the EKV to 
with K = 3 and the speed of the the ICBM for PG with K = 3. 
ICBM. 


The performance of the pursuit guidance can be increased by increasing the guid- 
ance coefficient, which might yield a smaller miss distance and intercept time. However, 
increasing the guidance coefficient will also increase the total commanded acceleration 
and cause an increase in the acceleration oscillation, making the EK V unstable. 


Cc. SIMULATION RESULTS FOR PROPORTIONAL NAVIGATION 
GUIDANCE 


In this section, the simulation results for PNG are analyzed. The first set of results 
are recorded when the acceleration command is calculated using the closing velocity. The 
second set of results are recorded when the acceleration command is calculated using the 
EKV velocity. 

1. Simulation Results for the Closing Velocity Approach 

The simulation is initialized with the same initial positions as the pursuit guidance 
in order to be able to make a fair comparison. The guidance coefficient N is selected as 5 
because both missiles have high acceleration forces acting on them. The acceleration lim- 
iter is set to 500 m/s? to prevent applying forces that the EKV cannot hold. The velocity 
vector of the EKV at the beginning of the interception is not towards the ICBM, which 
results in an unfavorable closing velocity. Hence, we picked an initial value for the clos- 


ing velocity as 6 km/s until the PNG steers the EKV towards the ICBM. 
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The three-dimensional interception geometry for PNG with N= 5 is shown in 
Figure 48. The EKV flies nearly directly to the collision point and hits the target. The 
EKV travels a curved trajectory at the beginning, which can be considered undesirable. In 
addition, the trajectory in endgame phase is also slightly curved due to the high accelera- 


tion of the ICBM. 


Space Launch Point 





Figure 48. Intercept geometry for PNG with N = 5. 


Figure 49 shows the miss distance for PNG with N = 5 as a function of scenario 
time. Figure 50 shows the altitudes of the EKV and ICBM. The minimum miss distance 
of this interception is 0.0182 m, which can be read from Figure 49. The minimum miss 
distance is clearly improved by using the PNG. The EKYV altitude is always above the 
ICBM altitude, which gives an advantage to the EKV, as shown in Figure 50. Staying 
above the ICBM altitude provides the EKV with a kinematical advantage over the ICBM, 
enabling the EKV to react to the ICBM’s evasive maneuvers. The interception time is 
read from Figure 49 as 3.25 minutes, which is within our intercept time window. The alti- 
tude of the EKV increases at the beginning, as shown in Figure 50, because of the initial 
angular velocity of the orbit, but the PNG starts steering the EK V downwards after ~0.7 


minutes. Note that this time was lower for the pursuit guidance. 
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Figure 49. Miss distance for PNG with Figure 50. Altitude of the EKV and the 
N=5. ICBM for PNG with N = 5. 


The minimum miss distance and the intercept time improved, but we have to 
question if we are reaching these results effectively. Hence, we also need to examine the 
total command effort used for the intercept. Figure 51 shows the instantaneous accelera- 
tion command of EKV as a function of scenario time under PNG with N = 5. Figure 52 
shows the cumulative acceleration. The instantaneous acceleration command reaches 169 
m/s? at the beginning of the intercept in order to put the EKV on the collision path, as 
shown in Figure 51. This acceleration decreases gradually as the EKV trajectory coin- 
cides with the collision path. The instantaneous acceleration presents a discontinuity at ¢ 
= 57s since it switches from the initial velocity value to the actual closing velocity to cal- 
culate the guidance command. The acceleration starts increasing at this point as a result 
of the change in the prediction of the collision path. The stage transitions can also be ob- 
served from the acceleration command of the EKV, as shown in Figure 51. Three minutes 
into the scenario, the guidance acceleration drops rapidly because the ICBM stops boost- 
ing and accelerating. The total acceleration command applied on the EKV can be read 
from Figure 52 as 21980 m/s*, which is considerably lower compared to the pursuit 
guidance. The oscillations observed in the instantaneous acceleration are a result of the 


autopilot. 
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Figure 51. Instantaneous acceleration ap- Figure 52. Cumulative acceleration applied 
plied on the EKV for PNG with on the EKV for PNG with N= 5. 
N=5. 


The velocities of both the EKV and the ICBM are shown in Figure 53. The veloc- 
ity of the EKV increases for 10s as a result of the single stage boosting. After the boost- 
phase, the velocity decreases slightly due to gravitational forces, until the EK V turns 
downwards, where the gravitational forces start acting in its favor. The velocity keeps in- 
creasing until the hit because the EKV always stays above the ICBM altitude pointing 


downwards. 


Figure 54 shows the closing velocity of the EKV to the ICBM as a function of 
scenario time for PNG with N = 5. The closing velocity of the EKV to the ICBM is small 
at the beginning of the intercept because the initial velocity vector of the EKV is not 
pointing to the ICBM. The PNG steers the missile towards the ICBM, which results in a 
rapid increase in the closing velocity, as shown in Figure 54. The closing velocity de- 
creases slightly as the ICBM velocity increases while boosting. The closing velocity stays 


constant after the third minute of the interception because the ICBM boost-phase ends. 
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Figure 53. Speeds of the EKV and the Figure 54. Closing velocity of the EKV to 
ICBM for PNG with N = 5. the ICBM for PNG with N = 5. 


The PNG improves the intercept time with a significantly smaller amount of com- 
mand effort, but it does not provide an optimum guidance at the beginning when we want 
the EKV to turn into the collision path as soon as possible. We will address this problem 
in the following sections. 

2. Simulation Results for Missile Velocity Approach 

The initial states of the EKV and the ICBM are the same as before for this simula- 
tion. We only alter the guidance calculation to use EKV velocity. For the missile velocity 
approach, the seeker is not required to measure the closing velocity. The guidance gain is 


set to 6 due to high acceleration of the EK V and the ICBM. 


The three-dimensional intercept geometry for PNG with N'=6 is shown in 
Figure 55. The initial performance is improved when compared to the closing velocity 
approach. The EKV turns towards the ICBM faster, and the curvature at endgame phase 


is decreased, which in turn decreases the intercept time. 
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Ge 
Figure 55. Intercept geometry for PNG with N'=6. 


Figure 56 shows the miss distance for PNG with N’ = 6 as a function of scenario 
time. Figure 57 shows the altitudes of the EK V and the ICBM. The minimum miss dis- 
tance for this simulation is 0.00095 m, which is a very successful result for a hit-to-kill 
intercept. The intercept time is decreased to 172 s and is within the 3.25 minutes window, 
which makes a 23.4 s improvement when compared to the closing velocity approach. The 
altitude of the EKV still increases before it compensates for the unfavorable initial veloc- 
ity vector, but this loss is also improved with the EKV velocity approach of PNG. The in- 
tercept takes place at an altitude of 250 km as shown in Figure 57. 
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Figure 56. Miss distance for PNG with — Figure 57. Altitudes for PNG with N’=6. 
N'=6. 
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Figure 58 shows the instantaneous acceleration applied to the EKV under PNG 
with N’ = 6. Figure 59 shows the cumulative acceleration. The initial acceleration in- 
creases to 325 m/s’, which puts the EKV on the collision path more effectively. The 
ICBM stage changes can also be observed from the EKV instantaneous acceleration 
chart. The PNG exerts more command effort on the EKV at the beginning of the inter- 
cept, which puts the EKV on the collision path faster, prevents extensive future guidance 
effort and decreases the intercept time. The total command effort is read from Figure 59 
as 21410 m/s, which is a 580 m/s’, improvement when compared to the PNG with the 


closing velocity approach. 
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Figure 58. Instantaneous acceleration for Figure 59. Total acceleration command for 
PNG with N’=6. PNG with N'=6. 


Figure 60 shows the speed of the EK V under PNG with N’ = 6 and the speed of 
the ICBM as a function of scenario time. Figure 61 shows the closing velocity. The ve- 
locity vector of the EKV is shown in Figure 60 and is very similar to the velocity in PNG 
with the closing velocity approach. The closing velocity increases faster when compared 
to the PNG with the closing velocity approach because the EK V turns more rapidly to- 
wards the ICBM. The maximum closing velocity is ~10.5 km/s for both PNG calcula- 
tions, but the EKV velocity approach reaches this velocity in 1 minute, where the closing 


velocity approach reaches this velocity in 1.7 minutes. 
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Figure 60. Velocities for PNG with Figure 61. Closing velocity for PNG with 
N'=6. N'=6. 


We have shown that the PNG using the EKV velocity approach performs better 
than the closing velocity approach. However, when we look at guidance unit’s perform- 
ance at the beginning of the intercept, we see that it still can be improved. The EKV is 
launched at an unfavorable attitude. By putting the maximum available acceleration 
command effort at the beginning of the scenario, the interception can be improved. Note 
that the current scenario is the most favorable we can achieve. We will face more unfa- 
vorable launch conditions when we have to launch the EKV at a later orbit time. We will 
address this problem in the next section. 

D. SIMULATION RESULTS FOR THE HYBRID GUIDANCE 

The proportional navigation guidance performs in the near optimal region for the 
window where the EKV is pointing toward the ICBM. The only non-optimal region is 
where the EK is in an unfavorable attitude due to initial orbital velocity. We need to 
steer the EKV towards the ICBM as soon as possible after the launch. To accomplish this 
goal, we will introduce the bang-bang guidance (BBG) method at the beginning of the in- 
tercept. BBG will be used until the ICBM enters the seeker’s gimble limit. The gimble 
limit of the EKV seeker is assumed to be 4.38 steradians. The guidance unit switches to 
the PNG after this point. This Hybrid Guidance (HG) method will provide the missile 


with a quick correction of the initial unfavorable attitude. 


The hybrid guidance simulation is run for the same conditions as the PNG, and 


the most significant improvement is achieved in the total command effort. The minimum 
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miss distance is 0.02854 m, the total intercept time is 170 s and the total command effort 
is 20,090 m/s? for the HG simulation. Note that the intercept time is improved by 2 s and 
the total command effort is improved by 1320 m/s? when compared to PNG with the 


missile velocity approach. 


Since we have reached the best guidance method that we can achieve for the 
given scenario, we will run this simulation for several EK V initial conditions and differ- 
ent ICBM launch points. In addition, there will be a launch delay due to target detection, 


launch decision and human interface delays as given by 
TS Ty + Tyo + Ty (4.43) 


where 1, is launch delay, 7, is detection time, 7, is the decision time and 1, is the human 


interface delay time. This is implemented in the simulation by delaying the EK V launch 


by 7, =30s. 


The three launch locations of interest are North Korea, China and Iran respec- 
tively. The initial conditions of the ICBM are defined by their initial state vectors. The 
EKV travels on a circular orbit with an altitude of 1000 km, an inclination angle of 
43.52° and a right-ascension angle of 15.28’. The initial state vector of the EKV on this 
orbit is defined by the mean anomaly @, which is measured from the ascending node. 
We will examine the guidance method for four different launch points within the launch 
range on the orbit. 

a) North Korea Case 
The initial state space vector of the ICBM is 


Xp -3029966.67 
Vp 3737980.67 


z, | | 4186406.53 
X,30)=| - |= (4.44) 
kp -279.69 
Dp 20.81 
Z 479.46 | cop 
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The minimum @ for an ICBM originating from North Korea is 99.97’ for a successful 


intercept within the desired intercept time window. The maximum q@ for the same case is 


114.72’. These angle values are determined after running the simulation several times. 


Table 8. Orbital launch points for the North Korea case. 





Position A B Cc D 











o 99.97° | 104.89° | 109.81° | 114.72° 























Figure 62 shows the four different intercept geometries against a North 





Korean ICBM for HG EKV guidance. The coverage angle is measured from point A to 
point D as 14.75’. A launch between these points on the orbit guarantees a successful hit- 
to-kill intercept before the ICBM releases its RVs. The intercept geometries of four dif- 
ferent launch points are shown in Figure 62. Any EKV carrier between point A and point 
D is capable of killing the North Korean ICBM. Before the EKV carrier leaves this re- 
gion from point D, another EKV carrier should enter from point A in order to assure the 
intercept of any ICBM launched from North Korea. Hence, we need a minimum of 


[360/14.75 | =25 EKV carriers in the orbit. 
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Figure 62. Intercept geometries against North Korean ICBM. 
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Figure 63 shows the miss distance of four different EK Vs under HG as a 
function of time. Figure 64 shows the altitudes of the EK Vs and the ICBM. The mini- 
mum miss distance for the EKV launched from point A is 0.01834 m. The EKV hits the 
ICBM 195 s after the ICBM launch, as shown in Figure 64. The significant improvement 
in the guidance method at the beginning of the intercept is observed in Figure 64, where 
the surplus altitude increase of the EK V is minimal compared to PNG results. The mini- 
mum miss distances and intercept times for points B, C and D are within the limits of the 
required constraints. Simulation results are listed in Table 9. Note that the time axes in 


the following figures start at the ICBM launch. 
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Figure 63. Miss distances with HG for Figure 64. Altitudes of EK Vs and North 
North Korean ICBM defense. Korean ICBM with HG. 





Figure 65 shows the acceleration command exerted on four EK Vs under 
HG as a function of scenario time. Figure 66 shows the cumulative accelerations. At the 
beginning of the intercept, the EK V uses BBG (as shown in Figure 65) until the ICBM is 
within the gimble limit of the EKV, which steers the EKV to the collision path as soon as 
possible. The linear increase region in the cumulative lateral acceleration is observed 
where the BBG is applied. The BBG is applied where the EKV is not pointing to the 
ICBM. Hence, the shorter the BBG applied, the better the Launch point. The preceding 
statement claims that point A is the most favorable launch point for an EKV against a 


North Korean ICBM. 
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Figure 65. Lateral acceleration exerted by Figure 66. Cumulative acceleration exerted 


HG on the EKVs launched from. by HG on the EKVs launched from 
points A, B, C and D against North points A, B, C and D against North 
Korean ICBM. Korean ICBM. 


Figure 67 shows the speeds of four EK Vs under HG as a function of time 
and speed of the ICBM. Figure 68 shows the closing velocities. Since the EKVs are 
launched with a 30-s delay, the ICBM has already attained a velocity of 555 km/s, which 
can be calculated using (4.44). The EKV speed increases rapidly for launch point C be- 
cause the collision path is perpendicular to the earth’s surface in this case. This case also 
introduces a head-on intercept since the EKV is literally right above the ICBM launch 
point. However, point D does not provide the EKV with an advantageous launch direc- 
tion. The HG exerts the maximum acceleration command on the EKV for 73 s to steer the 
EKV towards the ICBM, which causes the high total acceleration command observed in 
Figure 66. In addition, the closing velocity reaches its maximum value for launch point D 


because of the head-on intercept geometry. 
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Figure 67. Velocity magnitudes of the Figure 68. Closing velocity of the EKVs to 
EKYVs launched from points A, B, the ICBM for launch points A, B, C 


C, D and the ICBM with HG. and D with HG. 


The minimum miss distance, intercept time, intercept altitude and total 
command acceleration values are listed in Table 9; the intercept time is measured starting 
from the ICBM launch. The minimum miss distances are small enough for a hit-to-kill in- 
tercept. All the intercept times are before the ICBM delivers its RVs, and EKV is 


launched between point A and point D with a maximum launch delay of 30 s. 


Table 9. Simulation results for North Korean ICBM intercept with HG. 
































Point A Point B Point C Point D 
Miss Distance (m) 0.01834 0.01228 0.0839 0.1644 
Intercept Time (s) 195 145.68 138.06 194.16 
Intercept Altitude (km) 326.8 177.8 159.1 324.1 
Total Acceleration (m/s? ) 20110 21500 37330 26670 

















The LOS angular velocities for different launch points are shown in Figure 
69. The LOS angular velocity is proportional to the command acceleration where the 


PNG is applied. This proportionality proves the correct implementation of the PNG law. 
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Figure 69. LOS angular velocity magnitudes for points A, B, C and D with HG. 


b) China Case 
The initial state vector of the ICBM is 


X,(30)= 


T 
yr 


2r 
xP 


Jr 


_|3751990.63 


-985516.77 

5062656.41 

4.45 

-153.30 445) 
147.52 
533.47 


ECEF 


The minimum @ for an ICBM originating from China is 78.82° for a suc- 


cessful intercept within 3.25 minutes after ICBM launch. The maximum @ for the same 


case is 84.25”. 


Table 10. Orbital launch points for the China case. 





Position A 


B 


Cc D 











a 78.82” 














80.64" | 82.44° | 84.25° 














The coverage angle from point A to point D is 5.43’. A launch between 


these points on the orbit guarantees a successful hit-to-kill intercept before the Chinese 
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ICBM delivers its reentry vehicles (RVs). The intercept geometries of four different 
launch points are shown in Figure 70. Any EKV carrier that is present between point A 
and point D is capable of killing the Chinese ICBM. In this case, the minimum required 


number of EKV carriers is 67. 


Space Launch Points 





Figure 70. Intercept geometries against Chinese ICBM. 


Figure 71 shows the miss distance for four EK Vs under HG as a function 
of time. Figure 72 shows the altitudes of the EK Vs and the ICBM. The minimum miss 
distance for the EKV launched from point A is 0.08915 m. The EKV hits the ICBM 190 s 
after the ICBM is launched, as shown in Figure 71. The minimum miss distances and in- 
tercept times for points B, C and D are within the limits of the required constraints. The 
EKV hits the ICBM at an altitude of 287 km, as shown in Figure 72. Simulation results 
are listed in Table 11. Note that the time axes of the following figures start at the ICBM 


launch. 
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Figure 71. Miss distances with HG for Figure 72. Altitudes of EX Vs and Chinese 
Chinese ICBM defense. ICBM with HG. 


Figure 73 shows the acceleration command applied on the EK Vs under 
HG as a function of scenario time. Figure 74 shows the cumulative accelerations. At the 
beginning of the intercept, the guidance unit exerts the maximum acceleration on the 
EKVs (as shown in Figure 73) until the ICBM is within the gimble limit of the EKV 
seekers, which steers the EKVs to the collision path as soon as possible. Point A is the 
most favorable launch point for an EKV against a Chinese ICBM. The stage transitions 


of the ICBM are also observed in instantaneous lateral acceleration command. 
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Figure 73. Lateral acceleration exerted by Figure 74. Cumulative acceleration exerted 


HG on the EKVs launched from by HG on the EKVs launched from 
points A, B, C and D against Chi- points A, B, C and D against Chi- 
nese ICBM. nese ICBM. 
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The speeds of the EKV and the ICBM are shown in Figure 75. Since the 
EKV is launched with a 30-s delay, the ICBM has already attained 574 km/s of velocity, 
which can be calculated using (4.45). The EKV speed increases rapidly for point D 
launch, because the collision path is perpendicular to the earth’s surface in this case. The 
HG exerts the maximum acceleration command on the EKV for 76 s to steer the EKV 
towards the ICBM, which causes the high total acceleration command observed in Figure 
74. In addition, the closing velocity reaches its maximum value for launch point D be- 
cause of the head-on intercept geometry. The closing velocity of the EK V to the ICBM is 
shown in Figure 76. The closing velocity for point D launch reaches 22 km/s, which is 
approximately the sum of the EKV speed and the ICBM speed. This shows that the EKV 


is on a head-on intercept geometry. 








Speed in kms 











ae 1 15 2 25 3 35 
Time in minates 


Tice in minutes 


Figure 75. Velocity magnitudes of the Figure 76. Closing velocity of the EKVs to 
EKYVs launched from points A, B, the ICBM for launch points A, B, C 
C, D and the ICBM with HG. and D with HG. 


The simulation results for the Chinese ICBM interception from space 
launch points A, B, C and D are shown in Table 11; the intercept time is measured from 
the ICBM launch. The minimum miss distances that are listed in Table 11 are all smaller 
than the radius of the ICBM, which assures a hit-to-kill intercept. The intercept times 


show that the kills occur before the ICBM delivers its RVs. 
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Table 11. Simulation results for Chinese ICBM intercept with HG. 












































Point A Point B Point C Point D 

Miss Distance (m) 0.08915 0.09718 0.2521 0.2710 
Intercept Time (s) 189.96 171.42 153.42 124.26 
Intercept Altitude (km) 287.4 230.7 183.7 119.5 
Total Acceleration (ms”) 19780 19560 20010 26630 
Figure 77 shows the LOS angle rate as a function of scenario time. The 


HG guidance decreases the LOS angular rate during the intercept, as shown in Figure 77. 
If the LOS angular velocity was zero during the entire intercept we could say that the 
PNG was implemented perfectly. However, such perfect implementation is impossible to 
achieve for this scenario because the ICBM is accelerating, which is introduced to the 
guidance as a maneuvering target. Nevertheless, HG guidance performs extremely well 


and hits the ICBM with a very small miss distance. 
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Figure 77. LOS angular velocity magnitudes for points A, B, C and D with HG. 
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The initial state vector of the ICBM is 


Xp 3510071.07 
yp 4332137.25 
Bi 3095896.55 


X,(30) =|" |= (4.46) 
x kp 204.53 

Ve 43.94 

z 540.68 | ccnp 


The minimum @ for an ICBM originated from Iran is 38.35’ for a suc- 
cessful intercept within the desired intercept time window. The maximum @ for the same 


case is 47.57’. 


Table 12. Orbital launch points for the Iran case. 





Position A B Cc D 











a 40.01° | 41.92° | 43.44° | 45.73° 


























The coverage angle from point A to point D is 5.72’. A launch between 
these points on the orbit guarantees a successful hit-to-kill intercept before the ICBM de- 
livers its reentry vehicles (RVs). The intercept geometries from four different space 
launch points are shown in Figure 78. Any EKV carrier between point A and point D is 
capable of killing the Iranian ICBM. In this case, the minimum required number of EKV 


carriers is 63. 
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Figure 78. Intercept geometries against Iranian ICBM. 


Figure 79 shows the miss distances of four EKVs under HG as a function 


of scenario time. Figure 80 shows the altitudes of the EKVs and the ICBM. The mini- 
mum miss distance for the EKV launched from point A is 0.03729 m. The EKV hits the 
ICBM 190.5 s after the ICBM is launched, as shown in Figure 79. The EKV hits the 


ICBM at an altitude of 195.4 km as shown in Figure 80. Simulation results are listed in 


Table 13. Note that the time axes of the following figures start at the ICBM launch. 
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Figure 79. Miss distances with HG for Iran 
ICBM defense. 
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Figure 80. Altitudes of EKVs and Iran 
ICBM with HG. 


Figure 81 shows the instantaneous acceleration command applied to the 
EKVs by the HG as a function of scenario time. Figure 82 shows the cumulative accel- 
erations. The guidance unit applies BBG on the EKV (as shown in Figure 81) until the 
ICBM is within the gimble limit of the EKV, which steers the EKV to the collision path 
as soon as possible. Point D holds an unfavorable launch direction for the EKV, which 
requires it to sustain BBG for 94 s. Point D launch exerts the maximum total acceleration 
on the EKV as shown in Figure 82. Point A holds the most favorable launch direction for 
EKV against an Iranian ICBM. The stage transitions of the ICBM are also observed in in- 


stantaneous lateral acceleration command. 
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Figure 81. Lateral acceleration exerted by _‘ Figure 82. Total cumulative acceleration 


HG on the EKVs launched from exerted by HG on the EKVs 
points A, B, C and D against Iran launched from points A, B, C and D 
ICBM. against Iran ICBM. 


Figure 83 shows the speeds of four EK Vs under HG and the ICBM as a 
function of time. Figure 84 shows the closing velocities. The HG exerts the maximum ac- 
celeration command on the EKV where the ICBM is out of the gimble limit of the seeker, 
as shown in Figure 83. The EKV moves towards the earth’s surface in point D launch, 
which can be attributed to higher velocity. In addition, the closing velocity reaches its 
maximum value for launch point D as shown in Figure 84. For point D launch, the inter- 


cept geometry is a near head-on approach, which causes this high closing velocity. 
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Figure 83. Speeds of the EKVs launched Figure 84. Closing velocities of the EKVs 
from points A, B, C, D and the to the ICBM for launch points A, B, 
ICBM with HG. Cand D with HG. 


Simulation results are given in Table 13, where the time is measured start- 


ing from the ICBM launch. Note that there is a 30-s launch delay for the EK V launch. 


Table 13. Simulation results for Iranian ICBM intercept with HG. 


























Point A Point B Point C Point D 
Miss Distance (m) 0.03729 0.02465 0.02386 0.008616 
Intercept Time (s) 190.5 176.7 167.4 148.7 
Intercept Altitude (km) 195.4 168.6 152 123.4 
Total Acceleration (m/s” ) 27990 29310 20380 41410 




















The LOS angular rates are shown in Figure 85 and are proportional to the 


commanded acceleration where PNG is applied. 
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Figure 85. LOS angular velocity magnitudes for points A, B, C and D with HG. 


E. SUMMARY 

After running the simulation for North Korea, China and Iran, the performances 
of different guidance rules were analyzed, and a hybrid guidance is utilized for best per- 
formance. We showed that the hybrid guidance method performs extremely well for 
space-based ICBM defense for the scenarios of interest. We also showed that this guid- 


ance achieves a hit-to-kill intercept before the ICBM delivers its RVs. 
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VI. CONCLUDING REMARKS AND FUTURE 
RECOMMENDATIONS 


A. SUMMARY OF THE WORK 
This thesis developed a model for space-based, boost-phase interception of 
ICBMs from three example launch sites, namely North Korea, China and Iran. The ICBM 
was mathematically modeled and by taking the earth’s rotation and the atmospheric drag 
into account. The required orbit for a successful intercept was derived for the EKV carri- 
ers. Finally, a guidance algorithm was designed to achieve the interception for a hit-to- 


kill mission. 


The ICBM was modeled in three-dimensional space in the ECEF coordinate sys- 
tem. The initial launch parameters were calculated for a San Francisco, California attack 
by using the Lambert guidance. The precise parameters were generated by running a 


simulation of the developed model several times around the initial parameters. 


The orbit was derived to cover the three launch points of interest uniformly from 
space. The payload capability of the space launch vehicle and the lifetime of the system 
were also considered. The performance of the orbit was predicted using the maximum 
range of the EKV and the desired intercept features in the window of opportunity. 

The pursuit guidance (PG), proportional navigation guidance (PNG) and bang- 
bang guidance (BBG) are described and implemented in a three-dimensional model for 
the given scenario. The guidance rules were analyzed based on the simulation runs and 
their performances compared. The BBG and PNG were combined to form a hybrid algo- 
rithm for the best performance. 

B. SIGNIFICANT RESULTS 

Significant results obtained in this thesis work include the trajectory and the 

launch parameters of the example ICBMs, the required orbit for a space-based, boost- 


phase intercept, and an investigation of different guidance rules for the given scenario. 


The mathematical model we developed for the trajectory of the ICBM takes the 
earth’s rotation and the atmospheric drag into account along with the gravitational forces 


and the thrust acting on the ICBM. Since the ICBM travels in space during its mid- 
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course, the earth’s rotation plays a significant role on its impact point. For example, the 
impact point of the ICBM launched from North Korea with azimuth launch angle 39.6° 
and elevation launch angle 85.81° is San Francisco, California. However, if we take 
earth’s rotation out of the calculations, the impact point calculated by the simulation with 
the same launch parameters is N47° W 115.5°. The distance between San Francisco and 
this location is approximately 780 miles, which is the error if we assume a non-rotating 


earth. 


The orbit we derived covers the three example launch sites. The analysis made in 
Chapter III showed that the bounding limitation for the orbit altitude was the payload ca- 
pability of the space launch vehicle within the given scenario. The selected orbit has an 
inclination angle of i = 43.52°, a right-ascension angle of Q =15.28° and an altitude of 
ra= 1000 km. After running the intercept simulation it is found that the North Korea de- 
fense requires 25 EKV carriers in orbit, the China defense requires 67 EKV carriers, and 
the Iran defense requires 63 EKV carriers. In the interest of providing a common defense 
system, we have to consider the highest number, so the total required EKV carriers in the 


orbit is 67. 


Pursuit guidance (PG), proportional navigation guidance (PNG) and the bang- 
bang guidance (BBG) were implemented in a three-dimensional intercept model. The 
performance of these guidance rules were analyzed by running the SIMULINK® simula- 
tion for the North Korea case. From this analysis, the best performance was achieved by 
combining the BBG and the PNG guidance algorithms, which formed the hybrid guid- 
ance (HG). The HG is tested for three example launch scenarios, and four different space 
launch points were chosen for each launch scenario. The simulation results showed that 
the HG met the requirements of the space-based, boost-phase intercept as shown in Table 


9, Table 11 and Table 13. 


The most significant problem in the space-based interception is the initial launch 
direction of the EKV. In most cases, it is not possible to launch the EKV towards the 
ICBM due to angular velocity of the orbit. The initial launch direction of the EKV is cor- 
rected by using the BBG until the EKV pointed towards the ICBM. This method turned 
the EKV towards the ICBM as soon as possible by applying the maximum applicable 
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command acceleration and decreased the total command effort for the entire intercept. 
After ICBM entered the field of view of the EKV seeker, the PNG is used to steer the 
EKV, which minimized the command effort within the required intercept time. 
Cc. RECOMMENDATIONS FOR FUTURE WORK 

In this study, we assumed the earth as a perfect sphere and the ICBM as a point 
mass. These assumptions lead to inaccuracies in the ICBM trajectory calculation. In a fu- 
ture study, the oblate shape of the earth may be introduced to ICBM dynamics, and a full 


aerodynamic ICBM model may be considered for more accurate trajectory prediction. 


In this work, the EKV seeker is assumed to track the ICBM at all times and with- 
out any error. However, the EKV seeker has a limited field of view, so it cannot track the 
target at all times, especially at the beginning of the intercept. In a future study, the track- 
ing of the ICBM may be supported by the sensor network as illustrated in Figure 1 and 
investigated in previous studies [3] and [6]. The EKV seeker may also be remodeled by 


introducing the errors in tracking. 


In this thesis, we used the BBG and PNG as the selected guidance rules. Al- 
though these guidance rules perform satisfactorily, an optimum guidance algorithm may 


be investigated in a future study as proposed in [40]. 
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APPENDIX A. CODE FLOW CHART 


This appendix contains the flow chart of the MATLAB‘ code for the ICBM tra- 
jectory prediction and the SIMULINK“ model for ICBM intercept with different EKV 


guidance algorithms. 


START 
Main3d.m 
User decides 
launch location 
North lee, 




















China 
¥ ¥ 
Define N.Korea Define China Define Iran 
ICBM ICBM ICBM 
Parameters Parameters Parameters 
¥ 





Run Trajectory.m for 
selected launch 
location and angles 



























Calculate and 
save ICBM 
motion data 
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Read the 
given 
input data 


Assign initial 
values for looping 


Decide the stage with remaining 


propellant mass and reduce mass 
for given time step 


Define transition matrix F. 





Calculate next step state 
space vector using F 
X =X+FX 





Store state 
space 
vectors ina 
matrix 








Return 
calculated 
data to 
Main3d.m 
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START 
SimulationInit.m 


Calculate EKV location 
using OrbitDet.m 


Calculate EK 
state using InitMissile.m 


User defines ICBM 
launch location and 
launch dela’ 


Calculate ICBM initial state 
using 
TrajectoryForSimulnit.m 















Store Initial 
parameters to 
workspace 
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Simulate ICBM 
motion using /CBM 
Dynamics subsystem 


Track ICBM using 
Seeker subsystem 


Generate guidance 
command using 
Guidance subsystem 







or miss 





Steer the EKV with 
guidance command using 
EKV Dynamics subsystem 


| 
Store data in 
workspace 


Plot results using 
PlotGuidanceResults.m 
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Read remaining 
fuel and next 
state space vector 
from the 
integrators 


Decide the stage with remaining 
propellant mass 


Define transition matrix F, 


Calculate change in state space 
vector and the mass and send to 
integrator of Seeker subsystem 





Loop until EKV 
hits or misses 


Xdot = FX 
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Read ICBM and 
EKV data from 
ICBM Dynamics and 
EKV Dynamics 
sunsystems 


Calculate angles and pass to derivative 
tool of Seeker subsystem for deriving 
the angle rates 







Pass the 
track data to 
Guidance 
subsystem 
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Define coordinate 
conversion matrices 





Calculate guidance 
acceleration command in 
missile body coordinates 


Convert guidance 
command to ECEF 
coordinate system 


Pass guidance 
command to 
EKV Dynamics 
subsystem 
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Read remaining 
fuel and next 
state space vector 
from the 
integrators 









Define transition matrix F. 


Calculate change in state space 
vector and the mass 
Xdot = F X 






Loop until 

EKY hits or 

misses Add guidance acceleration 
command to state vector rate 









Pass the total 
state vector 
change to 

integrators 
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